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Rejection of Outliers* 


F. J. ANSCOMBE 
Princeton University and University of Chicago 
with collaboration by Irwin Guttman 
McGill University 


If one reading is a long way from the rest in a series of replicate determinations, or if 
in a least-squares analysis one reading is found to have a much greater residual than 
the others, there is temptation to reject it as spurious. Numerous criteria for the rejec- 
tion of outliers have been proposed and discussed during the past 100 years. They seem 
always to have been regarded as something like significance tests, and attention has 
been focussed on rejection rates. It is suggested that rejection rules are not significance 
tests but insurance policies, and attention would be better focussed on error variance. 
A detailed study is made of the effect of routine application of rejection criteria to rep- 
licate determinations of a single value. Determinations in triplicate and quadruplicate 
are especially considered. Complex patterns of observations are also considered, espe- 
cially factorial arrangements with high symmetry, and there is a study of the correla- 
tions between residuals. Attention is focussed mainly on rejection rules appropriate 
when the population variance is known, but some consideration is also given to Stu- 
dentized rules. 


1. Spurious OBSERVATIONS 


Variability or dispersion in a set of observations can arise from several different 
sources. Suppose that, for some reason we need not go into, it is desired to 
investigate the height (stature) of persons employed at a particular place. 
Three sources of variability in the readings are: 

(i) Inherent variability, the variability of stature that would be observed 
in the population even if all measurements were perfectly accurate. This vari- 
ability cannot be reduced without changing the population itself, the object 
of study. If we are interested in the mean stature of the population, we may 
refer to the variability as “error,” since it gives rise to estimation error. But 
the name is misleading. 

(ii) Measurement error, the error in using the measuring instruments. If all 
readings are made to the nearest centimeter, the measurement error should 
not exceed half a centimeter, but in fact it sometimes does. One may count 
as & measurement error any arithmetical mistake in reducing the original note- 
book entries to the form in which they are quoted as observations. 


* Research carried out in part at Princeton University and in part at the Department of 
Statistics, University of Chicago, under sponsorship of the Logistics and Mathematical Statis- 
tics Branch, Office of Naval Research. Reproduction in whole or in part is permitted for any 
purpose of the United States Government. 
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(iii) Execution error, a name we may use to include any discrepancy between 
what we intend to do and what is actually done, other than error in the use 
of measuring instruments. To include in the sample of measurements the height 
of some person not belonging to the population, to measure something other 
than height, or to select a biased sample—these would be classed as errors 
of execution. 

No observations are absolutely trustworthy. In no field of observation can 
we entirely rule out the possibility that an observation is vitiated by a large 
measurement or execution error. If a reading is found to lie a very long way 
from its fellows in a series of replicate observations, there must be a suspicion 
that the deviation is caused by a blunder or gross error of some kind. Several 
possible reasons why a reading might be grossly wrong can usually be thought 
of without difficulty. In such cases, the reading will be checked or repeated 
if that is possible. If not, it may be rejected as spurious because of its big resi- 
dual, even though there is no other known reason for suspecting it. In suffi- 
ciently extreme cases, no one hesitates about such rejections. The question 
is, where should the line be drawn? One sufficiently erroneous reading can 
wreck the whole of a statistical analysis, however many observations there are. 

Statements to the above effect have appeared many times in the literature. 
Even writers who have expressed total disapproval of the rejection of outliers 
may be found to insert a parenthetical remark, ‘‘except for obviously incorrect 
readings.” Reasonable argument can occur, not over the permissibility of ever 
rejecting outliers, but over the completeness of the rejection, how far the fact 
of a rejection can be forgotten—as well as over the question of where pvecisely 
to draw the line. If we could be sure that an outlier was caused by a large measure- 
ment or execution error which could not be rectified (and if we had no interest 
in studying such errors for their own sake), we should be justified in entirely 
discarding the observation and all memory of it. The act of observation would 
have failed; there would be nothing to report. Such an observation could justly 
be described as spurious. If, on the other hand, we could be sure that an outlier 
was caused, not by any large error, but by some peculiarity (nonnormality) 
of the inherent variability of the population under study, then it might still 
make good sense to discard the observation from a statistical analysis based 
on the method of least squares, but the observation should not be forgotten. 
A correct statistical summary of the observations would include a report about 
outliers as well as an analysis of the remaining observations (cf. Kruskal, 1960). 

If a rather pronounced degree of nonnormality is expected, it may be advisable 
to transform the observations, or to work with medians instead-of means, or 
some other such device.* A modified least-squares method, with weights depend- 


* The suggestion about medians in this context goes back to Edgeworth (1887). It is not 
always satisfactory in principle to make use of any such device. If the subject of study is the 
output of an agricultural or industrial process, where most of the variability is inherent, the 
mean yield is inescapably what ought to be measured, not the median yield or mean logarithmic 
yield (for instance). Perhaps it should be explicitly stated that in this paper we are concerned 
only with observations whose primary purpose is the estimation of location parameters of some 
sort. Observations made primarily for the measurement of dispersion are not considered. 
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ing on the residuals, has been suggested by Jeffreys, with particular reference 
to seismology and astronomy. Although such a technique is possible as a guard 
against spurious observations, it will not be considered here. We envisage a 
situation where the intrinsic variability and ordinary measurement errors result 
in a nearly enough normal paitern of variation for simple averages or equal-weighted 
least-squares estimates to be attractive, bul where it is feared that occasionally an 
observation will be affected by a gross error. 

An observation with an abnormally large residual will be referred to as an 
outlier. Other terms in English are “wild,” “straggler’”’, “sport” and ‘‘maverick’’; 
one may also speak of a “‘discordant,” “anomalous” or “‘aberrant’’ observation. 
Spurious will mean: affected by an abnormally large measurement or execution 
error. Usually it is not possible to check directly whether an observation is 
spurious; it is a matter of conjecture. 

The object of this paper is to study the effect of unthinking routine use of 
a specified rejection rule on the estimation of means (or linear combinations 
of means). There are circumstances where the problem of spurious observations 
can appropriately be met by an impartial rejection procedure, which is thought 
of as part of the experimental technique. It is not suggested that the problem 
should always be met this way, but the results of the present study should 
be helpful (though not decisive) to anyone who wishes to judge individual 
outliers “on their merits,”’ intelligently. 

The gist of the paper can be had by reading as far as Section 3, skimming 
Sections 4 and 5, and then turning to Section 10. Section 7 on complex designs 


is to some extent a digression, and may have interest apart from the subject 
of outliers. 


2. SKETCH OF THE History OF REJECTION OF OUTLIERS 


The subject of rejection rules has been held to be important and interesting 
in many sciences—astronomy, geodesy, chemistry, physics, ballistics. Nowa- 
days it may be thought of as lying within the broader subject of “data proces- 
sing,” an essential though perhaps too little considered branch of statistical 
analysis. It has been discussed in countless books on the combination of obser- 
vations by least squares, over a period of nearly a hundred years, as well as 
in many contemporary books on statistics. 

One of the first references to the rejection of outliers seems to have been 
a remark by the leader of the German school of astronomers, Bessel, in a geo- 
detic work published in 1838, to the effect that he had never rejected an obser- 
vation merely because of its large residual; all completed observations, with 
equal weight, ought to be allowed to contribute to the result. ‘‘We have believed 
that only through strict observance of this rule could we remove arbitrariness 
from our results.” 

The first attempt at a rejection criterion based on some sort of probability 
reasoning was that of Peirce (1852). Peirce’s argument was reproduced by 
Chauvenet (1863), who then gave a similar rule based on a simpler argument. 
A heated discussion was provoked by these rules, and by some other suggestions 


put forward for disposing of outliers. In due course the topic became standard 
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in books dealing with least squares; lengthy historical surveys were given by 
Czuber (1891) and Wellisch (1909). See also Rider (1933). 

The next significant step after Peirce and Chauvenet was apparently made 
by Wright (1884). After an excellent general discussion of the problems raised 
by the occurrence of outliers among astronomical readings, he suggests that 
the best rule for a computer to follow who is not the observer is to reject any 
observation whose residual exceeds in magnitude five times the probable erro: 
(ie. 3.37 times the standard deviation). The reason given for this is that if 
the Gaussian law of error is truly satisfied, only about one observation in a 
thousand will be rejected, ‘‘and therefore little damage will be done in any 
case.” 

From 1925 onwards statisticians have paid much attention to the subject. 
Student (1927) gave an interesting account of analytic shemical determinations, 
and proposed the use of a range criterion. Thompson (1935) may be said to 
have “Studentized” Wright’s rule. A good survey of this literature has been 
given by Grubbs (1950). A recent development is the Bliss-Cochran-Tukey 
rule (1956), developed for use with several of the bioassays in the U. 8. Pharm- 
acopeia. 

Rejection of outliers seems to have been a peculiarly American topic. At 
any rate, Peirce, Chauvenet and Wright were Americans, and most of the 
recent developments have occurred in the U. S. The subject arose in the context 
of least squares, or, as current jargon has it, Model I analysis of variance; and 
it is curious that (apparently) it was not considered by Gauss. 

Chauvenet, after he had reproduced Peirce’s argument, introduced his own 
simplified version of the rule by some remarks beginning: ‘“The above investi- 
gation of the criterion involves some principles, derived from the theory of 
probabilities, which may seem obscure to those not familiar with that branch 
of science.” The same may be said of Chauvenet and. of everyone else; and 
familiarity with that branch of science does not remove the obscurity. It is 
easy now to laugh at Peirce and Chauvenet, and at other eminent nineteenth 
century figures such as Glaisher and Bertrand, as they wallow in the probability 
quagmire. But have today’s statisticians been so much clearer headed? All 
published proposals for rejection criteria, based on any kind of mathematical 
reasoning, from Peirce’s onwards, have an unexplained starting point or objec- 
tive, presented as though it were the only obvious one and in fact utterly obscure. 

All writers seem to have thought of the problem as something in the nature 
of significance testing.* They have varied in the significance level judged ap- 
propriate for action, but all have regarded significance levels (or something 
of the sort) as relevant beyond discussion. Peirce and Chauvenet thought in 
terms of something like a 50% significance level per set of data. If applied to 
a simple sample of n observations drawn from a homogeneous normal source, 
either criterion rejects as often as not, very roughly, so the rejection rate is 


* Jeffreys puts this the other way round and suggests that the general use of tail-area prob- 
abilities by statisticians for significance testing goes back to Chauvenet’s criterion (Jeffreys, 
1939, p. 316; 1948, p. 357). 
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something like 1 per 2n observations. Wright’s rule rejects at something like 
the rate of 1 per 1000 observations, almost independently of n. Most modern 
statisticians, bemused by 5%, give rules having rejection rates of about 1 per 


20n observations. No one has explained why this should be so. No one seems 
to have asked. 


3. THe NATURE OF REJECTION RULES 


Rejection rules are not significance tests. In a study of whether, and how 
often, clearly spurious observations occur in a certain field, significance tests 
may be appropriate. But when a chemist doing routine analyses, or a surveyor 
making a triangulation, makes routine use of a rejection rule, he is not studying 
whether spurious readings occur (he may already be convinced they do some- 
times), but guarding himself from their adverse effect. The same is often true 
of a statistician when he performs least-squares analyses, especially if he has 
a single very large bulk of data or many similar smaller sets of data to analyze. 
The statistician may feel little temptation to exercise personal judgment in 
“processing” the data, and prefer to follow a rigid rule for treating outliers, 
whose long-term effect he knows. 

A rejection rule is like a householder’s fire insurance policy. Three questions 
to be considered in choosing a policy are 

(1) What is the premium? 

(2) How much protection does the policy give in the event of fire? 

(3) How much danger really is there of a fire? 

Item (3) corresponds to the study of whether spurious readings occur in 
fact—a study that is hardly possible unless plenty of readings are available. 
The householder, satisfied that fires do occur, does.not bother much about (3), 
provided the premium seems moderate and the protection good. In what cur- 
rency can we express the premium charged and the protection afforded by a 
rejection rule? That depends on the purpose of the observations; an answer 
can be given as soon as a suitable loss function is specified. 

Often estimation errors are only one among several types of error affecting 
the whole investigation or process. For example, the error made by the manu- 
facturer in assaying a drug is not the only reason why the patient receives an 
incorrect dose. The doctor judges only roughly how much of the drug the patient 
needs, the prescription is filled by the pharmacist with less than infinite pre- 
cision, and the quantity consumed by the patient may bear only a faint re- 
semblance to the instruction of one teaspoonful every four hours. So long as 
the assay errors are on the whole small, there would be no interest in very 
detailed information about their distribution; the bias and variance will suffice. 
In research experiments, the immediate object may be to measure responses 
to certain treatments, but the ultimate object is usually to throw light on some 
broad class of phenomena. Estimation errors are not the only impediment to 
clear perception. If the investigation is in the end abortive, estimation errors 
will most likely not be solely responsible. Again, we do not need to know more 
about the estimation errors than their variance (and bias, if any). 

So variance will be considered here, although in principle any other measure 





128 F. J. ANSCOMBE 


of expected loss could be used. The premium payable may then be taken to 
be the percentage increase in the variance of estimation errors due to using 
the rejection rule, when in fact all the observations come from a homogeneous 
normal source; the protection given is the reduction in variance (or mean 
squared error) when spurious readings are present. Rejection rates are of no 
more than incidental interest. 

With this approach, any suggested rejection rules can be investigated and 
compared. As in all previous published work on outliers, it will be assumed 
here that the following conditions are satisfied (they need not be in practice). 

(i) Whatever circumstance causes an observation to be spurious is not expected 
to affect neighboring observations; all observations are supposed independent 
in this respect. It is also supposed that spuriousness is uncorrelated with the 
reading that would have been obtained had the observation been made without 
abnormal error. : 

(ii) Computation costs can be ignored. If that were not so, the “premium”’ 
would have to include the extra computation cost resulting from using the 
rule. 

(iii) No prior knowledge concerning the means or regression coefficients that 
are to be estimated from the data is incorporated in the rejection rule, which 
is therefore “impartial.’”” For example, suppose the observations consist of just 
two replicate readings, the population mean is to be estimated and the popula- 
tion variance is known. Any impartial rejection rule must lead to rejection if 
the difference between the readings exceeds some critical value, and then since 
there is no way of saying which is the better observation both readings must 
be rejected (or one after flipping a coin, perhaps). In practice, the observer 
may consider one of the readings reasonable and the other unlikely, and so 
he will retain the first and reject the other. He will then be using prior information 
about the population mean. 

This restriction to an impartial rejection rule has an important effect on 
the character of our results. We return to consider it further in Section 7. 

It is natural to consider rejection criteria based on the magnitude of the 
residual. Such a type of criterion seems likely to work better than other relatively 
easy criteria based on order statistics. In the least-squares analysis of complex 
patterns of observations, the computation of residuals is, or should be, a standard 
procedure, and a rejection criterion based on residuals is therefore particularly 
convenient. To begin with, we consider the simplest possible case, appropriate 
to typical chemical analyses, where only one mean is to be estimated and the 
population variance can be supposed known. We then pass on to complex 
patterns, and also consider ‘‘Studentized”’ criteria for use when the population 
variance is unknown. 

Notation. Throughout, the observations will be denoted by y, , yz , --* 5 Yn; 
there being n observations in all. The residuals (before rejection of any outliers) 
will be denoted by z, , 22, °** , 2, . v will always denote the number of residual 
degrees of freedom, i.e. v is the rank of the matrix transforming (y;) to (z,). 


4. FORMULATION OF THE PROBLEM FOR A SIMPLE SAMPLE 


We are given observations y; , ¥2 , °** » Yn (XW = 3). It is hoped they are a 
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random sample from a normal population N(u, o”), where o is known and yu 
is to be estimated. But possibly one or more of the y; are spurious, coming 
from a different source, and ought to be rejected. We consider the effect of 
applying a rejection rule routinely to samples of fixed size n. Let 


2=%y4-79, ng = Dy (@¢@ = 1,2, +++ ,n). 


If y; is omitted, the average of the remaining observations is 


2» yi/v = 9 — z:/v, (4.1) 


where y = n — 1. More generally, if several observations are omitted, say 
Yr» Y2»°** » Yr, the average of the rest is 


G9—@tat::: +2,)/nm—7). (4.2) 


Let M be the serial number of the observation having the greatest residual, 
so that 


lzu| > |z.| forall i+ M. (4.3) 


(We suppose that the observations are recorded to sufficient decimal places 
for no two residuals to be equal in magnitude.) 

We propose to reject any observation whose residual is excessively large. 
The following type of rule is unsatisfactory: 


RULE 0. For given C, reject every observation y; such that 
lz;| > Co . 
Estimate u by the mean of the retained observations. 


The reason is that a single outlier, if it outlies sufficiently, can cause all the 


\z)’s to exceed Co, and the whole sample would then be rejected. The following 
rule is more cautious: 


RULE 1. For given C, reject yy if \zu| > Co; otherwise no rejections. Estimate 
u by the mean of the retained observations, thus 


R=7 af zar| < Co, 
=%-—zu/v tf |eu| > Co. 


Under this rule, not more than one observation can be rejected. The most 
frequent values for n in chemical analysis are 3 and 4, and then Rule 1 is prob- 
ably the best that can be suggested. If more than one observation out of a 
very small sample appeared to be spurious, the observer would most likely 
wish to scrap them all. For large samples, however, the possibility of multiple 
selective rejections needs to be considered, and the following rule is suggested: 


RULE 2. Apply Rule 1. If an observation is rejected, consider the remaining 
observations as a sample of size n — 1 and apply Rule 1 again; and so on. Esti- 
mate wp by the mean of the retained observations. 
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It would be possible for the values of C to differ in the successive applications 
of Rule 1, but there is no obvious advantage in this, and, in so far as Rule 2 
is considered below, C will be supposed constant. In fact, it is difficult to study 
Rule 2 exactly, apart from Monte Carlo computation. Rule 1 is easier, and 
sometimes has almost the same effect as Rule 2—namely, when there is not 
more than one spurious observation present in the sample, C is not very small 
and n is not very large. 


5. THrory ror Rute | (Simpie SAMPLE) 


To study the long-run effect of application of Rule 1 to samples of size n, 
we consider the distribution of 2, when the y; are interpreted as chance vari- 
ables rather than as a particular realization of chance variables. The z; and M 
become chance variables too. 

If there are no spurious observations, the joint distribution of (2; , 22, -** , 2n) 
is independent of the distribution of 7, so 7 and zy are independent. Each z; 
has variance vo"/n; (n/v)* z;/o has the standard normal distribution N(0, 1). 
Let a random variable 7 be defined as the following function of zy : 


T=0 if leu| < Co, 
—(n/v)*(eu/o) if \e| > Co. 
Then a = 9 + oT /(nv)', 9 and T are independent, &(@) = u and 
ae ( sf) 
var (f) = = 1+ ds (5.2) 


The rejection rate (proportion of observations rejected in the long run) is 


(5.1) 


= ch flew| > Co} =ch (T £0}. (5.3) 


If an observation is spurious, it is convenient to think of it as being equal 
to an observation from N(u, o”) plus an extra independent error or bias. If 
there is one spurious observation among a sample of n, and if the bias is large 
enough, the observation is almost certain to be rejected under Rule 1 or Rule 2. 
If there are two spurious observations, and their biases are large and the dif- 
ference of their biases is large, it is almost certain that both will be rejected 
under Rule 2. 

Suppose a particular observation, y, say, is from N(u + be, o”), while 
Y1 » Y2 > *** » Ya-1 are independent from N(u, o*). Then @ is distributed 
N(u + ba/n, o”/n), independently of (2, , 2. , --* , 2). Under Rule 1, 


a = {9 — bo/n} + {oT /(rw)* + bo/n}, (5.4) 
where T' is defined as before by (5.1). We have 


sa — u)? = “(1 re * g(r + bVoh), (5.5) 


and the rejection rate is still given by (5.3). It is to be observed that the distri- 
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bution of zy , and therefore that of 7, depends on 6 and is different from what 
it was before. In the absence of any rejection rule, (5.5) would be replaced by 


eg - == (142), (5.6) 


In principle, (5.2) and (5.5) can be evaluated by quadrature, though if n 
is higher than 4 some skill in classical n-dimensional geometry is called for. 
The joint distribution of the (z,;) is spherical normal in a »-dimensional flat 
space. If there are no spurious observations the distribution is centered at the 
origin; if one (or more than one) observation has an added bias the center 
is away from the origin. No observation is rejected provided that 


—Co K.2i°< Co 


for all 7, so that the point (z,) lies inside a region R bounded by n pairs of parallel 
(v — 1)-flats equidistant from the origin. For n = 3, v = 2, R is the interior 
of a regular hexagon; for n = 4, v = 3, R is the interior of a regular octahedron. 
If (z;) lies outside R, which observation is rejected, and what sign its residual 
has, are determined by which one of the 2n faces of F is intersected by the line 
segment joining (z;) to the origin. To evaluate (5.2) or (5.5) an integration 
must be performed over the exterior of R. For (5.2) there is great symmetry, 
and &(T’) can be expressed, after one integration along the radius vector, 
as 2n times an integral over any one face of R. The faces are regular simplexes, 
and the integrand is a function only of distance from the center of the simplex. 
Thus for n = 3, v = 2, each face is a line segment, and we obtain 


T 


ae 6 uve Hessen 2 1 ) dt 
am)=" fe C+ 


itei+e 


For n = 4, v = 3, each face is an equilateral triangle and the integrand is a 
function of distance from its center, rather more cumbersome than the above. 
It is far more difficult to deal with the case of one spurious observation, since 
the center of the distribution no longer coincides with the center of R, and 
much of the symmetry is lost. Only n = 3 has been considered. Numerical 
results obtained are given below in Section 10. For further progress in the 
accurate evaluation of (5.2) and (5.5) it will no doubt be necessary to resort 
to Monte Carlo techniques. 


6. APPROXIMATE FoRMULAS 


Fortunately, it seems that we do not have to wait for accurate calculations 
to obtain results good enough for practical use. Asymptotic approximations 
can be found that are easy to use and ought to be tolerably reliable, as a guide 
to action. Just how reliable they are can only be shown convincingly by com- 
paring them with correct values. The comparisons in Section 10 are highly 
encouraging, as far as they go. 

For the case of no spurious observations, the proportion of observations 
such that the residual exceeds Co in magnitude, which we shall denote by a, 
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is given by 


a = 26(—(n/r)'C), 


where 


a) = [ oy du, $0) = 6°/V2m. 


If 1. denotes the two-tailed a-point of the standard normal distribution, we 
have , 


t. = (n/r)*C. (6.1) 


a would be exactly the rejection rate if we used Rule 0. If C is very large, a 
is very small. Moreover, if C is large, residuals exceeding Co in magnitude will 
tend to occur singly; it will be relatively unusual for two or more observations 
in the same sample of n to have residuals exceeding Co in magnitude. This 
is easily proved by noting that the joint distribution of any pair of residuals 
is normal with correlation coefficient different from +1. By taking C large 
enough, we can make the conditional chance that |z;| > Co, given that |z;| > Co, 
as small as we please, and hence the conditional chance that any one of |z;| > Co 
for j = 1, 2, --- , n, with j ¥ 7, given that |z;| > Co, as small as we please. 

It follows that, asymptotically as C — ©, we can ignore the possibility that 
more than one residual in the sample will exceed Co in magnitude, and we 
obtain 


&(T’?) ~n / f(t) dt 


tl>te 
= n{2t.o(ta) + a}, (6.2) 


and a@ is asymptotically the rejection rate under Rule 1. 

If C is not so large as effectively to rule out the possibility that more than 
one of the |z,;| exceeds Co, but nevertheless we substitute (6.2) into (5.2), it 
may be expected that the result will approximate the effect of Rule 2 rather 
than Rule 1, especially if n is fairly large. We obtain this result for the simple 
sample if we suppose the rejected residuals to be independent and replace the 
divisor n — r in (4.2) by n — 1. 

Table 1 shows 2¢,¢(¢,) + a and a as functions of t, . The entries were readily 
obtained from the National Bureau of Standards table of probability functions. 
For comparison with the results of Section 11 it may be noted that. (6.2) can 
be expressed alternatively in terms of the tail area of a x’ distribution with 3 
degrees of freedom. 

Suppose now that there is just one spurious observation, y, , drawn from 
N(u + be, o”), and let us suppose that b is fairly large. Then it is very probable 
that M = n. Let us assume this to be so, for the moment. Then it is easy to 
see that 7 + (v/n)*b has the following distribution: outside the interval 


= ((v/n)'b — (n/v)*C, (/n)'b + (/)*C) 
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Tasie 1 
Table of 2t,¢(t,) + a, with a shown in brackets 


0 8 
03826 03630 03442 .03263 03092 


(00373) (.00350) (.00328) (.00308) (.00288) 


.02929 .02773 02625 02483 02348 
(.00270) (.00253) (.00237) (.00221) (.00207) 


.02219 02096 01980 .01869 01763 
(.00194) (.00181) (.00169) (.00158) (.00147) 


.01663 .01568 .01478 .01392 01311 
(.00137) (.00128) (.00120) (.00111) (.00104) 


01234 .01161 01091 01026 00964 
(.00097) (.00090) (.00084) (.00078) (.00072) 
00905 00850 .00798 .00748 .00701 
(.00067) (.00063) (.00058) (.00054) (.00050) 
00657 00616 00577 00540 00505 
(.00047) (.00043) (.00040) (.00037) (.00034) 


00472 00442 00413 00385 .00360 
(.00032) (.00029) (.00027) (.00025) (.00023) 


the distribution is that of a standard normal variable N(0, 1); inside I the 
chance is concentrated at the mid-point (v/n)'b. We therefore have 


&(T' + (v/n)*b)? = 1 + / (vb?/n — #)g(t) dt. 


For b positive and large, we can replace the upper end-point of I by ~, obtaining 


&(T + (v/n)*b)? = 1 + (b’/n — 1)®(—2) — x¢(z), (6.3) 
where 


x = (v/n)*b — (n/r)*C. 


If (6.3) is substituted into the right-hand side of (5.5), the result will be a 
lower bound for &(@ — «)*, because we have ignored the possibility that M may 
be different from n. We may hope, however, that the bound will be close enough 


to the true value to indicate how large b needs to be for the rejection rule to 
give good protection. 


7. CompLex PATTERNS oF DaTA 


Let us turn now to complex patterns of data. We suppose the observations 
Y. » Y2» *** » Yn (if none is spurious) to be drawn from independent normal 
distributions having common variance o and means that are given linear 
functions of some unknown parameters. At first we shall suppose o” known. 
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When all the unknown parameters have been estimated by least squares, the 
residuals {z,} can be calculated. 

We observe first of all that in general the residuals do not all have the same 
chance distribution. This is illustrated by simple linear regression: y; is drawn 
from N(u + 62; , o”), where the x; are predetermined, and let us say 2,2; = 0. 
Then 


2=4%-9- (By sxj/2jx4)a; ’ 


and we find 


var (z;) = (2 a i Jet. 


n Da; 


Thus 0 < var (z,) < (n — 1)o’/n, both bounds being attainable for some i 
and some set of {z,}. If the {z;} are equal-spaced, we have 


(n — 4)o"/n < var (z;) < (n — l)o’/n, 


and if n is not very small we shall not go far wrong if we suppose all the resi- 
duals to have the same variance, say vhe average value over 7, which is 
(n — 2)o”/n, or vo"/n. All the residuals have-exactly this variance if (and only if) 
there are just two different z-values, with n/2 observations at each. 

_In order to keep our considerations as simple as possible, we shall restrict 
attention to cases where (in the absence of spurious readings) all residuals 
have exactly the same variance, namely vo"/n. All ordinary factorial designs, 
where the different levels of each factor are replicated equally often, have 
this property, and also Latin squares and balanced incomplete block designs. 
A type cf design to be excluded is a factorial design where the levels of one 
of the factors are not equally replicated, as in an agricultural experiment on 
insecticides, where it would usually be advisable to replicate the “control” 
treatment of no insecticide more heavily than any one of the alternative in- 
secticide treatments under test. It is to be expected that results obtained under 
the assumption of equal variances will be approximately correct if the resi- 
duals have only approximately equal variances, but we shall not examine this 
question. 

We observe next that in general not every pair of residuals has the same 
correlation. For some designs, some of the correlation coefficients are +1. An 
example is the 3 X 3 Latin square, where the residuals are equal in sets of 3, 
according to the letters of the orthogonal square, and so the correlation coeffi- 
cients between the 36 possible pairs of the 9 residuals are 9 of them equal to 1 
and 27 of them equal to —}. Suppose we are given that ¢ = 1 and the lay-out 
and yields are as follows, where 0, 1, 2 denote the three levels of a treatment: 


O n=129 () »=59 (2) %=63 
(2) y% = 6.0 (0) ys = 5.7 (1) ys = 6.4 
(1) y= 6.0 (2) ys = 63 (0) y% = 4.9 
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The residuals are 


4 >= 2.04 = —1.22 23 —0.82 
@ = —1.22 = -082 2= 2.04 
2, = —0.82 = 2.04 2 = —1.22 


Since the standard deviation of residuals is (v/n)'s = 0.47, the values 2.04 
may be judged excessively large and indicative of a spurious reading. The 
problem is, which reading? If we reject any one of y, , Ys and ys and treat that 
observation as missing, the estimated value to be substituted is equal to the 
observed value minus n/vy (= 4.5) times the residual (2.04), i.e. we subtract 
9.2 from one of y; , ¥s , Ys , and proceed to estimate row, column and treatment 
effects in the usual way. The new residuals are 


4 = Cc 23 = —0.20 23 0.20 
a4 = —0.20 25 0.20 26 0 
2, >= 0 .20 23 = 0 ~a=- 0.20 


and these are satisfactorily small. But our estimates of the row, column and 
treatment effects will depend sharply on which one of the three questionable 
observations is rejected. If we insist on an impartial rejection rule, not using 
any prior knowledge concerning the unknown parameters to be estimated, we can- 
not have any preference, and to settle the deadlock we might consider rejecting 


all three observations. We should then scrap the whole experiment, because 
(in the absence of prior knowledge) too few readings are left for the unknown 
parameters to be uniquely estimated. 

In practice the experimenter always has some prior knowledge, which may 
be small compared with good observations having ¢ = 1 but is perhaps ap- 
preciable compared with a gross error, here estimated at 9.2. If he rejects y. 
or. ys , the fitted yield will be negative, which he may find disturbing. Apart 
from that, he very likely thinks that rather small row, column and treatment 
effects are more probable than large ones, and therefore the estimates obtained 
after rejecting y, are more plausible than those obtained after rejecting either 
Ys Or Ys . In fact, he may guess that just the initial digit of y, is spurious, and 
y, should read 3.9. 

One may conclude from this example that it is unwise to try to apply an 
- impartial routine rejection rule to data such that some pairs of residuals have 
correlation +1. It may be wise to go further and exclude designs such that 
some pairs of residuals have large correlations, say greater than } or # in magni- 
tude, because if a spurious reading occurs it is not unlikely that another reading 
will be rejected instead. Thus it is of interest to know, in respect of any proposed 
design, what correlations the residuals will have, and this is a matter not discus- 
sed in books on the design of experiments. 

Let the observations {y;} be represented by the » X 1 column vector y, 
let @ be a (n — v) X 1 vector of nonredundant unknown parameters and let A 
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be the n X (n — v) matrix of known coefficients, such that, in the absence of 
spurious readings, 


y has a spherical normal distribution with &(y) = A®é. (7.1) 


A is of rank n — », and V exists such that A’A = V™. The n X 1 vector of 
residuals z, after least-squares estimation of 6, is given by 


z = Qy, (7.2) 


where Q = I, — AVA’, I, being the n X n unit matrix. Q has the following 
properties: 


(i) Q isan X n symmetric idempotent matrix of rank v. It correctly trans- 
forms y to z, when 6 is estimated by least squares as if (7.1) were true, whether 
or not (7.1) is true in fact. It is invariant to transformations of the parameter 
set 6. 

(ii) The 7th row of Q shows how a gross error in y;, is distributed among 
the residuals (as we see by replacing y by a vector consisting of 1 in the 7th 
place and 0’s elsewhere). 

(iii) The variance matrix of the residuals, if (7.1) is true, is given by 


&(zz’) = Qo’. (7.3) 
Two practical methods of calculating Q are:- 


(i) Solve the least-squares equations in any convenient manner and express 
the {z;} in terms of the {y;}. For sufficiently symmetrical designs, this need 
be done for only one z, say 2, . 

(ii) If it is convenient (it is always possible) to choose ® so that V is a scalar 
multiple of I,_, , say V = vI,_, , then Q = I, — vAA’, and to find AA’ we 
calculate scalar products of pairs of rows of A. This method is easy for 2* fac- 
torial designs—for which, however, it is almost the same as (i). 

If all residuals have the same variance, all coefficients in the principal diagonal 
of Q are equal to v/n, and from the relation Q? = Q we easily deduce that 
the mean squared correlation coefficient between pairs of residuals is 


SF taken: 
, (n — 1) 


Thus the square root of this quantity is a lower bound to the magnitude of 
the largest correlation coefficient between any pair of residuals. Only rarely 
(i.e. for a small proportion of available designs) is this lower bound attained. 

If the observations have a simple two-way classification, in k rows and / 
columns, with n = kl, and the expectation of each observation y is a row con- 
stant plus a column constant, the possible correlations between residuals are 
easily seen to be 


(7.4) 


1 


1 
-? > 


(7.5) 


where «x = k — 1,4 = 1 — 1. If either k or 1 = 2, there are correlations of —1. 
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If the observations have a one-way classification, say k observations in each 
of l columns, so that the expectation of each observation is just a column con- 
stant, the possible correlations between residuals are 


oom oO; 


K 





(7.6) 


again we have correlations of —1 if k = 2. In particular, there will be cor- 
relations of —1 if the design consists of two replications, in either one block 
or two blocks, of some factorial design and if all interactions between the factors 
are estimated, as well as the general mean and main effects and (for two blocks) 
the block difference. 

Table 2 lists the correlations that occur among the residuals for some designs 
with fairly small n. Designs having correlations of +1 have been omitted, 
except for designs (2) and (14), which have been included for comparative 
purposes. The correlations for design (2) are an example of (7.6); those for 
design (3) an example of (7.5). In designs (1) and (2), three two-level factors 
are compared in eight observations. Both designs can be obtained by deleting 
four factors from a “saturated” 2’/16 design, as given by Plackett and Burman 
(1946); and it is interesting to note that two essentially different designs can 
be obtained in this way. Design (4) is obtained by deleting any six factors 
from the Plackett-Burman saturated fraction of 2'*. Designs (7) and (8) are 
half replicates of a 2° factorial. (7) is the obvious choice, and achieves equal- 
magnitude correlations. (8) illustrates the effect of another alias subgroup. 
The design ‘2°/2: ABC” has correlations differing from those of (8) in sign 
only. Designs (9), (11) and (14) were given by Brownlee, Kelly and Loraine 
(1948) and have often been quoted since, as the “optimum” fractions, such 
that main effects are not confounded with any two-factor interactions. With 
respect to correlations among residuals they are not optimvm, and are bettered 
by designs (10), (13) and (15), which, with (16) and (17), have alias subgroups 
permitting the lowest possible largest correlation. Design (12) is of interest, 
because it can be reinterpreted as two replications (in two blocks) of a 2° fac- 
torial, with all two-factor interactions estimated. If the three-factor interaction 
is also estimated, we have at once the correlation pattern of (7.5), with some 
correlations equal to —1. Designs (16) and (18) are essentially the same design; 
like (7), they achieve equal-magnitude correlations. Design (24) can be re- 
interpreted as a 3° factorial with all two-factor interactions estimated. Of all 
these twenty-four listed designs, (4) alone is unsymmetrical, in the sense that 
all the rows of Q do not contain exactly the same set of coefficients. 


8. ForRMULATION FOR CoMPLEX PATTERNS 


Restricting attention, then, to designs such that, in the absence of spurious 
observations, all residuals have the same variance and no two residuals have 
correlation coefficient equal, or very close, to +1, we wish to examine the 
effect of routine application of a rejection rule. Two questions to consider first 
are: what sort of rule, and effect on what? 

Let M be defined as before, by (4.3). Rules 0 and 1 of Section 4 can be adapted 
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TABLE 2 
Correlations between residuals for various designs 


Correlations 
between 
Design residuals 


(1) 2° ; +4, 0 

(2) 28/2 repeated: ABC ; —1,0 

(3) 3? 

(4) Plackett-Burman fraction of 2° 

(5) Balanced incomplete blocks 

(v= 4,k = 3) 

(6) 2¢ 

(7) 25/2: ABCDE 

(8) 25/2: ABCD 

(9) 28/4: ABCD, CDEF 
(10) 2/4: ABC, DEF 
(11) 27/8: ABCD, CDEF, ACEG 
(12) 27/8: BCE, CDF, BDG 
(13) 27/8: ABC, CDE, EFG 
(14) 28/16: ABCD, CDEF, ACEG, EFGH 
(15) 28/16: ABC, CDE, EFG, AGH 
(16) 29/32: ABC, CDE, EFG, AGH, BFI 
(17) 2°/64: ABE, ACF, ADG, BCH, 

BDI, CDJ 0.383 
(18) Latin square 6 0.333 2+} 
(19) Balanced incomplete blocks 
(v 2K 7, k= 3) 0.285 —4, i, =} 
(20) Latin square 0.212 —i,% 
(21) Graeco-Latin square 0.298 3, -} 
(22) 3 0.116 ° —t,y45 —-a 
(23) 34/3: ABCD 0.139 +4, 0 
(24) 39/729: ABD, AB*E, ACF, AC*G 
BCH, BC*I 0.302 —43,4, -} 


Most of the designs are factorials. Thus (10) is a quarter replicate of a factorial arrangement 
with 6 factors A, B, C, D, E, F, each at 2 levels; ABC and DEF are generators of the alias 
subgroup. (2) is « half replicate repeated in one block. (4) is given by Plackett and Burman 
(1946). In the balanced incomplete block designs (5) and (19), v is the number of varieties or 
treatment levels, & is the number of units in each block. It is to be understood that no interactions 
are estimated. Several designs can be reinterpreted as fewer-factor designs with two-factor 
interactions estimated, as mentioned in the text. For the incomplete block designs, interblock 
information is supposed not recovered. 


to the present case, by changing the second sentence to read: Estimate the un- 
known parameters from the retained observations by the method of least squares. 
Normally it will be appropriate to perform this estimation by the two-stage 
method known misleadingly as ‘‘analysis of covariance,” the rejected observa- 
tions being regarded as missing values to be estimated. 

Suppose a particular observation, y, say, is missing or (if present) is held 
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to be spurious. Let {z;} be the residuals when some arbitrary value for y, is 
used in the least squares estimation of the parameters. Then if y, is replaced 
by y¥: — 2, the residuals are replaced by {z$"}, where 


a) 
zi 


= 2% — 19» (8.1) 
{q:;} being the coefficients of Q. The value of » which minimizes 2,(z;'’)’ is 
§ = Veida/Bqh - 


Given our assumptions, we have 2.q3, = qi = v/nand 22g = 2ijQagdiiy; = 
Zihiys = 2, . Thus 


4 = (n/v)z, . (8.2) 


If all observations except y, are good, the correct least-squares estimates of 
the unknown parameters are found by substituting y, — 4 for y, in the expres- 
sions that would have been valid if all observations, including y, , had been 
good. The residuals are now given by (8.1) with (8.2) substituted, that is 


qa) a) . 
a =0, ee Fh Pah (i ~ 1), 


where p;, is the correlation between z; and z, ; and we obtain 
var (2?) = (1 — pis)vo"/n. (8.3) 


Thus Rule 1 (which is the only rule that we consider in detail) prescribes 
that if |z~| > Co, ys must be replaced by yy — nz,/v in the parameter estimates. 

How should our former Rule 2 be adapted? We have just seen that only 
for some rare designs are all the correlations p;; equal in magnitude. (8.3) shows 
that if one observation is rejected, the new residuals will have unequal variance, 
in general. Provided the correlations are fairly small, it may seem reasonable— 
it is certainly simplest—to take no account of the changes in variance in formu- 
lating the rule, which would run: Apply Rule 1. If an observation is rejected, 
compute revised residuals and apply Rule 1 again; and so on. Finally compute 
the least-squares estimates of the unknown parameters from the retained observations. 

Consider now how to measure the effect of routine application of one of these 
rules. We have argued in Section 3 above that it is appropriate to consider 
the sampling variances of estimates. Which estimates? In a complex design 
we are not necessarily interested in all the parameters that have to be esti- 
mated. The general mean and the block differences (if any) are usually of no 
more than secondary interest; that is, the experiment cannot be said to have 
been performed in order to estimate them. Sometimes the same is true even 
of the responses to some of the treatment factors. Let us therefore divide the 
parameters into two sets, those that are “‘of interest’’ and the remainder, and 
propose to study the variance matrix of the estimates of the parameters “of 
interest.”” (We shall impose a mild restriction on the choice of parameters for 
the “of interest” set.) Let the number of parameters “of interest’? be x. By 
an orthogonal transformation, the covariances of the parameter estimates can 
be made zero, and the product of the variances after transformation is equal 
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to the determinant of the variance matrix before transformation. A possible 
measure of the “premium” charged by a rejection rule is the percentage in- 
crease in the determinant of the variance matrix caused by using the rule, 
when in fact no observations are spurious. We shall see that, for Rule 1 with 
given C, this percentage increase is proportional to x. Rejecting an observation 
implies a loss of information concerning each of the parameters (in the absence 
of spurious observations), and in a sense it is true that the greater « is the greater 
is the loss. 

However, to aggregate the losses for each of the x parameters is not altogether 
reasonable as a measure of the premium exacted by the rule—as we see by 
considering the grouping of observations. It is commonly held that to group 
normally distributed observations at a grouping interval not exceeding o/2 
has a negligible effect on the estimation of means, and is therefore to be recom- 
mended for ease of computation, whenever ease matters. Such a grouping 
increases the residual variance by about 2%, and so the sampling variance 
of every estimated parameter by the same percentage. If the reader considers 
that the recommendation about grouping is equally sound however many 
parameters are to be estimated, he will presumably consider the following 
definition of the premium charged by a rejection rule to be reasonable: the 
premium is that proportional increase in o° which would inflate the determinant 
of the variance matrix of the estimated parameters of interest by as much as the 
rejection rule does, when no observations are spurious. 


9. THEoRY FoR RuLE 1 (CoMPpLEX PATTERNS) 


We are supposing that (in the absence of spurious readings and with no 
rejection rule applied) 


(i) all residuals have the same variance, and 

(ii) no two residuals have correlation coefficient equal or close to +1. 

We shall also suppose that 

(iii) if just the x parameters of interest are estimated, the remainder being 
set equal to 0, the resulting residuals all have equal variance. 


These conditions concern the design—specifically, they condition the linear 
space spanned by the columns of A—and they also concern the choice of the 
x parameters of interest out of the total of n — v parameters. Despite references 
to variance and correlations, they are independent of distribution assumptions. 
Condition (iii) implies, for example, that we are not allowed to pick out as the 
“of interest” set a one-degree-of-freedom component, e.g. the linear. component, 
from the main effect of a three-level factor. 

Let the parameters of interest be the first « components of 6. We wish to 
find the ratio of determinants of the variance matrices of the estimates 
(6, , 6, , --+ , 6.) calculated under the assumptions (a) that Rule 1 is applied, 
(b) that Rule 1 is not applied—or that Rule 1 is applied with C = ©. We observe 
first of all that Rule 1 does not depend on the parameterization. If a nonsingular 
lower-triangular transformation T is applied to @ and the problem is rephrased 
in terms of 6* = T6, no change is made in the residuals or in the operation of 
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the rejection rule, and T transforms the estimates of @ to the estimates of 6*, 
ie. 6* = T6. In particular, the leading x-rowed minor of T transforms 
(6, , 6. , «++ , 6,) to the corresponding estimates of interest (6% , 64 , --- , 6). 
It follows that the ratio of determinants that we seek is the same for the latter 
estimates as for the former. There is no loss of generality in supposing that 
the parameters 6 satisfy the conditions 


(iv) the parameters of interest are the first x components of 6, 
(v) V =L_,/n. 


(If (v) is not already satisfied, apply the lower-triangular transformation T 
to 6 defined by nT’T = V™.) 

From (i) and (v) it easily follows that the sum of squares of the coefficients 
in any row of A is equal to n — »; (v) also implies that the columns of A are 
orthogonal and the sum of squares of coefficients in any column is equal to n. 
If we now introduce condition (iii), it is easy to deduce that the sum of squares 
of the first « coefficients in any row of A is equal to x. In the absence of any 
rejection rule, the least-squares estimates of 6 are 


6 = A’y/n. 


Rule 1 has the effect of replacing yy by yu + (n/v)*cT, where T is defined 
by (5.1). Let us consider distributions conditionally on the value of M. We 
can apply an orthogonal transformation to the parameters of interest, so that 
the transformed A satisfies 


Gyo = Qy3 = *** = Qy, = 0, 


and then, because the transformation does not disturb the sum of squares of 
the first « coefficients in any row of A, we must have 


ay = Vr. 


(The sign is arbitrary; let us take the positive root.) The rejection rule does 
not now affect 6, , 6, , --- , 6, , but 6, is modified by the rejection rule by the 
addition of the term 

(x/nv)*oT 


If in fact no observations are spurious, this added term is independent of 
what it is added to (which has variance o”/n), and so we get the following result. 
We can transform the parameters of interest, by an orthogonal transformation 
depending on M, so that one new parameter has its variance changed to 


(o*/n){1 + («/r)8(7)}, 


while all the other components have unaltered variances of o”/n. The determi- 
nant of the variance matrix (which is invariant to orthogonal transformation 
of parameters) is thus 


(o*/n)‘{1 + («/r)8(T")}. (9.1) 
This result is independent of M and therefore true unconditionally. The factor 
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in curly brackets is the ratio that we are looking for, and is invariant to any 
nonsingular transformation of the parameters. The premium, according to the 
suggested definition at the end of Section 8, is 


{1 + «/)a(T*)}"* — 1, (9.2) 
which is less than, but may be closely equal to, 


(1/)8(T"), (9.3) 
unless x = 1, when these expressions are exactly equal. 

Now suppose that one observation, say y, , is spurious and has bias be. It 
may happen that M = n, or, if not, that the first « coefficients in the Mth row 
of A are the same as those in the nth row. Conditionally on this being so, our 
previous line of reasoning shows that the determinant of the matrix of mean 
squared estimation errors is changed from (9.1) to 


(o?/n)*{1 + (&/v)8(T + (v/n)*)*}. (9.4) 


This result does not hold unconditionally, but suffices for application of the 
results of Section 6. If no rejection rule were applied, (9.4) would be replaced by 


(o*/n)"{1 + (/n)b’}. (9.5) 


The arguments of Section 6, used to derive (6.2) and (6.3), are still valid 
in this context. We are now in a position to draw some conclusions. 


10. Numericat Resutts 


Both for the simple sample and for the fairly broad type of complex patterns 
just considered, we see, from (6.2) and (5.2) or (9.3), that the premium charged 
by Rule 1 can be roughly reckoned at 


(n/») {2t.0(te) + a}, . (10.1) 


where i, = (n/v)*C, as explained in Section 6. This should be tolerably accurate 
if C is so large that na is quite small. If C is not so large, there is some reason 
to expect that (10.1) approximates the premium charged by Rule 2 rather 
than Rule 1. 

How much premium we are willing to pay should depend on how greatly 
we fear spurious observations. But, as with domestic insurance, we shall probably 
not care very much provided the premium is small. Let us see what can be had 
for a premium of 2%, i.e. an effective increase of 2% in the residual variance. 
Setting the expression (10.1) equal to 0.02, we obtain the values of C and a 
given in the first part of Table 3. The lower part of the table refers to a premium 
of 1%. The value 1 for v/n is not attainable, but a simple sample of large size 
has v/n close to 1. The least value of v/n for the designs given in Table 2 was 
8/27, or about 0.3. Lower values for v/n are no doubt rare, except perhaps 
in some screening experiments. 

The protection given by Rule 1, when one observation is spurious with a 
large bias be, is given approximately by substituting (6.3) into (5.5) or (9.4) 
and comparing the result with (5.6) or (9.5), respectively. The best the rule 
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could do would be always to reject the spurious observation, and in that case 
we should have 


&(T + (v/n)*d)? = 1. 


This is the limit of (6.3) as b — ©. To give some idea of how large b must be 
for the rule to work well, values of b are tabulated for which, according to (6.3), 


&(T + (v/n)*b)? = 1.5. 


The chance that the spurious observation will escape rejection is, in the cases 
tabulated, of the order of 0.02. For somewhat smaller b the rejection rule will 
be beneficial, better than no rejection rule, but the mean squared estimation 
error will be substantially greater than if there had been no spurious observation. 
If b = (n/v)C, the chance that the spurious observation will be rejected is 
about 0.5. 

For given n and no spurious observations, as v is reduced the effect of a re- 
jection becomes more serious. Consequently the rejection rate a must be reduced 
to preserve a fixed premium, and (n/v)'C is increased. On the other hand, 
for a fixed bias b, (v/n)*b is reduced, and so the z in (6.3) is reduced twice over. 
Hence the striking increase in the b’s as one reads across Table 3. If we increase 
the premium a little, the b’s are not greatly reduced. 

Thus for an experiment of fixed size n, as the design is made more ingenious 
and v becomes smaller, a gross error in one of the readings spoils more estimated 
effects and becomes less detectable. Even if we can prevent the correlations 
from blowing up to +1, it takes a bigger error to be seen. 

The above findings indicate that, while the rejection rate a should depend 
on »/n, it should not otherwise depend upon n. To that extent they support 
Wright against most other authors, who have advocated fixed rejection rates 
per experiment. 

It remains to compare the approximate formulas with exact calculations. Com- 


TABLE 3 
Tabulation of rejection rules, based on approximate formulas 


v/n 1.0 0.8 0.6 0.4 


2% premium 
Cc 3.14 2.56 1.63 
a 0.00171 0.00094 0.00026 
b 5.1 5.8 6.9 8.7 12.8 


1% premium 
Cc 3.37 3.08 2.73 2.31 1.72 
a 0.00076 0.00058 0.00042 0.00026 0.00012 
b 5.4 6.1 7.2 9.1 13.3 


For given premium, values of C and a are given, and also that value of b such that a bias 
of size be in one observation inflates the mean squared estimation error only one and a half 
times as much as rejection of one observation at random when all observations are good. 
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TABLE 4 


Tabulation of rejection rules, based on accurate calculations 


Approximate 
Rejection Approximate rejection 
C (n/c?) var (f) rate (n/o*) var (f) rate 


Simple sample, n = 3, v = 2 

2.46003 1.04 0.002433 1.04241 0.002588 

2.66184 1.02 0.001065 1.02088 0.001114 

2.84623 1.01 0.000475 1.01032 0.000490 

3.01724 1.005 0.000214 1.00512 0.000220 

Simple sample, n = 4,» = 3 

2.57994 1.04 1.04134 0.002891 

2.79541 1.02 1.02043 0.001247 

2.99206 1.01 1.01014 0.000550 

3.17434 1.005 1.00504 0.000247 

The table refers to Rule 1 when there are no spurious observations, and shows values of C 
corresponding to given premiums of 4%, 2%, 1%, and 4%, together with the rejection rate. 


The last two columns show for comparison values obtained by the approximate formulas, 
using the C values quoted. 


[Results for n = 3, v = 2, one spurious observation, are to follow.] 


putations for the simple sample with n = 3, Rule 1, were begun by Mr. William 
W. Hardgrave, using a desk machine. Later, one of us (I. G.) used an IBM 650 
to check and extend the calculations. Some results are presented in Table 4. 


11. c UNKNOWN 


Clearly the above investigation leaves many questions unanswered. Further 
calculations are needed to determine the reliability of the approximate formulas. 
One would like to know the effect of routine application of a rejection rule to 
data such that the component of random variation was homogeneous but not 
normal—for example, having a distribution of the same shape as Student’s 
distribution with 7 degrees of freedom (Pearson’s Type VII with exponent 4), 
as suggested by Jeffreys, or having a “contaminated” distribution of the type 
considered by Tukey (1960). How does the rule function in the presence of 
Tukey’s (1949) removable nonadditivity, or of other systematic departures 
from standard assumptions? 

In view of statistical history during the past fifty years, the question likely 
to be asked first is: what happens if ¢ is unknown? It is not harder’in principle 
to investigate Studentized rejection rules than those already considered. Let 
s’ denote the best (quadratic unbiased) estimate of o° available. If derived 
from the given observations (y;) only, s” has v degrees of freedom. But in general 
we may suppose that there is also some prior information about o”, equivalent 
to a quadratic estimate having », degrees of freedom, so that s° has »v + 
degrees of freedom. The condition |zy| > Co in Rule 1 and in the definition 
of T,, (5.1), will be replaced by |zy| > Cs. 
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For the case of no spurious observations, we can obtain an approximate 
expression for the premium, corresponding to (10.1), as follows. Choosing any 
one of the residuals, z, say, for consideration, we can make a fixed orthogonal 
change of axes so that (z;) is transformed to (f, , {2 , --- , f, , 0, --- , 0), where 
where the ¢’s have independent N(0, o”) distributions, and ¢, = (n/v)*z, . It 
follows that the condition |z,| > Cs is equivalent to 


(HoH) _ a) > ot, (11.1) 


where x” stands for a x” variable with » + v) — 1 degrees of freedom, inde- 
pendent of ¢, . The chance a that this condition is satisfied, which is approxi- 
mately the rejection rate under Rule 1, is given by 


f= Vv a , rms 1] Pa gi7tvent) 
y vty —nC’/v ii 


where the right-hand side denotes the two-tailed a-point of Student’s distri- 
bution with v + ») — 1 degrees of freedom. This formula corresponds to (6.1) 
and extends a result given by Thompson. It can be expressed alternatively 
in terms of the incomplete beta function ratio tabulated by Karl Pearson, thus: 


v + a i } 
q = pe ae ae ? 2 ? 
where x = 1 — nC’/v(v + »). The premium »~*&(7”) is given approximately 
by (n/vo”) times the partial expectation of ¢{ , integration being confined to 
the region where (11.1) is satisfied. We obtain the result: 


1 gq) 8 7(Pt%=1, 3) 
; &(T”) 5 ls 3 5) (11.3) 
where z is as above. This may be expressed alternatively in the following rule 
for determining C, given the premium. Multiply the premium by »v/n, and 
find the corresponding upper percentage point of the variance ratio (F) distri- 
bution with 3 and v + ») — 1 degrees of freedom. Calling this F, we have 


nc? ee 
vy =-1 + BF — 1)/0 +») 


For example, suppose that v/n = 0.5, vy + v = 30 and the premium is to be 


0.02. We find that the upper 1% (= 0.02 X 0.5) point of F with 3 and 29 degrees 
of freedom is 4.54 and so 


(11.2) 


13.6 
2C* = a 9.59, C = 2.19, 


and from (11.2) we have a = 0.00092. If » + ») = 121, the premium and v/n 
remaining the same, we find similarly C = 2.33, a = 0.00079; and for vy + » = © 
(¢ known) we have C = 2.38, a = 0.00076. For a given premium, C depends 
quite sharply on v + vo , and to a lesser extent a also varies with »v + ». 

For such numerical calculations, Federighi’s table (1959) is useful for (11.2) 
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and Table 16 in the Biometrika Tables (1954) for (11.3). Ingenious interpolation 
is called for. 


It would be possible to investigate the response of the rejection rule to a 
spurious observation, as in Section 6, but I have not done this in general. The 
extreme case is the simple sample with n = 3, v = 2, » = 0. Rule 1 is now 
the same as Rule 0, because at most one observation in the sample can differ 
from the mean by more than C's (when C > 1); and it follows that formulas 
(11.2) and (11.3) give exact results for Rule 1. For a premium of 2% we find 
we must take C = 1.154638, and then a = 0.00667. The rejection rule can be 
expressed in the form: reject either extreme observation if its distance from the 
median observation exceeds 82 times the distance of the other extreme obser- 
vation from the median. If one of the three observations has a bias of magnitude 
bo, the other two being good, the chance that the rejection rule will caus? the 
exclusion of the bad observation is about 3 if b = 79; and b needs to be three 
or four times this size before rejection of the bad observation can be said to 
be virtually certain. For all practical purposes the rejection rule is utterly 
useless and absurd (cf. Lieblein, 1952). One may conjecture that a Studentized 
rejection rule will have low power whenever v + » is small, say less than 30. 
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TECHNOMETRICS May, 1960 


Locating Outliers in Factorial Experiments 


CuTHBERT DANIEL 
New York City 


Editor’s Note: The following paper by Mr. Cuthbert Daniel was presented at 
the Washington, D. C. Meetings of the American Statistical Association on December 
29, 1959, as part of a session on the subject of outliers. Mr. Daniel later became aware 
of an error (see footnote associated with Equation 6) that may invalidate some of 
the conclusions reached in the paper. Nevertheless, the Editor has decided to publish 
the paper so that a complete record of the meeting is available, and in the belief that 
the paper represents an important beginning in the study of patterns of residuals 
resulting from carefully plan experimental designs. 


1. InTRODUCTION 


An outlier in a factorial experiment is an observation whose value is not in 
the pattern of values produced by the rest of the data. The pattern has two 
aspects, and two kinds of disturbance of the pattern may appear. That observa- 
tion is suspect whose removal greatly simplifies the description of the rest of 
the data. That observation is also suspect that complicates the description of 
the error distribution. 

An outlier may be due to a single relatively large interaction in one cell; 
it may be an extreme manifestation of the random error system; or it may be 
a mistake. Unless they are quite large, unreplicated factorial experiments can 
not be counted on to supply positive evidence for deciding between these three 
alternative possibilities. Partially replicated factorial experiments permit an 
unbiased and independent error estimate, and so put the experimenter in a 
better position to spot bad values and to decide on their source. 

An advantage of full replication (I believe not mentioned by Fisher) is the 
security against outliers given by having at least two independent observations 
for each treatment combination. The empirical cumulative distribution of the 
ranges of the replicate pairs is often all that is necessary to detect a bad value 
and to decide that a mistake has been made (Mistakes are in my experience 
the commonest source of bad values. Such mistakes are usually not clerical 
but technical.) When the error distribution is normal, the ranges of pairs have 
the half-normal distribution and so can be plotted about a straight line on 
half-normal probability paper (1). The magnitude of the largest range or two 
can then be put in proper perspective. 

Bad values are, then, easily found in fully replicated experiments. A strong 
case for such replication can be made when there is prior evidence of frequent 
bad values. But this evidence is often not available and the expermienter may 
then decline to spend much of his time in replication. In my judgment the 
statistician should always urge at least partial duplication (2, 3), especially 
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in the smaller fractional replicates of the 2?~* form when (p — q) is 4 or 5. In 
larger experiments, and when the run-to-run variance is securely known, single 
replicates may be all that can be managed. This paper deals mainly with un- 
replicated experiments. 


2. ONE BIASED VALUE IN A Two-Way Layout 


The row and column classifications are here assumed discontinuous or quali- 
tative. The first panel of Table 1 shows a set of (synthetic) data in a 5 X 4 
layout. The second panel shows the “regression values’ computed from the 
row averages and the column deviations under the assumption of additivity. 
The third panel gives the residuals, that is, the differences between the original 
values of the first panel and the corresponding regression values of the second. 

The second row of panels in Table 1 shows the same set of data, handled 
in the same way, but with one value (third row, third column) decreased by 20. 
It is seen that the largest residual is in the cell (3, 3). However it is not 20 less 
than the former residual but only 12 less. 

The third row of this table shows the effects of this one change in the 5 X 4 
layout. The last row in Table 1 shows the effects on the residuals of a bias of 
size RC in a single cell (here chosen as (1, 1)) of an R X C layout. Here r is 
R — 1 and cis C — 1. It is obvious that the same pattern can be centered on 
any cell. 

A biased value will produce a pattern in the calculated discrepancies, d;; , 
that should show high correlation with the constants, c;; , given at the end of 
Table 1, when the latter are centered on the maximum discrepance, d,,,. . The 
numerator of the correlation coefficient, 


(1) Vesidis = (Lo 05; d;)/( Qo sy dis)” 
is, as always a contrast. We can re-write equation 1 
(2) a ee 6,,/\weRC)(SSD) 
where SSD is the sum of squares of the discrepances and 
a R Cc R Cc 
(3) 6, = - Dei d;; = Zz De Cis Yis = RC dusx 
i=1 j=1 t=1 j=1 


is the ‘“‘maverick contrast’’. The third member of equation 3 is useful for esti- 
mating the variance of the maverick contrast because the observations, y;; , 
are often safely assumed statistically independent. Thus Var (6,,) = o” >.>. ci; . 
The last member reveals that all the information about a single bad value is 
in the largest residual, d,,,. , after all. Using the identity 


R Cc 
(4) DX Lie; = reRC, 
‘ i 
the variance of the largest discrepance is seen to be rc/RC times the error variance. 


Substituting the values given by equations 3 and 4 in equation 1, and abbre- 
viating the mean square for all RC residuals by MSD, we get 
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(5) 2 we &.. a 
e1i-4is “" re SSD rc? MSD 


which shows how tr” varies with the size of the layout and with the observed 
discrepancies. This formula for t® shows the fractional reduction in the residual 
sum-of-squares, SSD, that is possible by revision of the response in the cell 
showing the largest residual. For the data of Table 1, as modified in the second 
set of panels, this fraction is 0.669. 


TABLE 1 
Effects of one biased value on residuals in an additive two-way layout 


“Data” Regression Residuals 


35 32 40 33 -1 -2 
29 29 36 29 0 -2 
25 29 40 28 3 
19 25 35 23 3 
22 20 29 22 —2 -2 


Avs. SB 2 & 
Devs. -—-4 -—3 6 


Biased Data Biased Regression Biased Residuals 


35 32 40 34 39 38 —2 
29 29 36 30 35 34 —1 
25 29 20 24 29 28 5 
19 25 35 24 29 28 1 
22 20 29 23 «28 = «27 —3 


26 27 32 
—-3 -2 3 


Bias in Data Bias of Regression Bias of Residuals 


0 —3 -1 -1 3 -1 
0 —3 -1 -1 3 -1 
—20 —8 4 4 -—12 4 
0 —3 -1 -1 3-1 
0 —3 -1 -1 3-1 


—4 
—3 


0 


*>con gw CO oooeoce 
Omo ecceco 
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If we had only the spoiled data, and did not know that 20 had been removed 
from the observation at (3, 3), we would have spotted its location from the 
large discrepance. We would then have estimated a simplifying value to enter 
in its place. The easiest way to do this is to revise the response by the amount 
—RC/rc times the discrepance d,,,, . For the numbers given the change would 
be plus 15. Thus the present 20 would be changed to 35. The new SSD is (1 — r’) 
times the cld SSD. For our case the new SSD is then (1 — 0.669)202 or 67. 
The new estimate of sigma, called s, , is (67/11)! or 2.47. The earlier estimate, 
based on the biased residuals would be (202/12)? or 4.10. 

Since the contrast 6,, is chosen after the data are taken, it seems natural to 
use a modification of Scheffé’s S-method (pages 68-72, reference 4) for judging 
the reality of the bias, and hence of the largest residual. It will be assumed in 
this section that a fair estimate of the true error variance can be obtained from 
the revised discrepances, with (rc — 1) degrees of freedom. This estimate, 
called s? below, can be computed from Equation (6)’. 


(6) 3 = (ssp - ae a2.) / (re — 1) (see footnote) 


Scheffé’s method permits us to find critical magnitudes of the maverick 
contrast, that is, values that will be judged just significant,'at some chosen 
confidence level. Using primes on 6, and on dmx to denote just-significant 
values, we place 


(7) On = Sé(bn); 

where 

(8) S = (reF,)', 

and F, is the (1 — a)% critical value of F, that is 

(9) Fy = Fo—e),re.re-1 3 

and where the estimated standard error of the maverick contrast is given by 
(10) 6(6,,) = 8,(reRC)*. 


Substituting the values given by equations 4, 8, 9 and 10 in equation (7) gives 
the following simple relationship between d/,,, and s, : 


(11) RC dyax = 8,(reF'o)*(reRC)*. 
Defining /,,, as the ratio of these two quantities as shown in equation 12, 


(12) bia = Anex/8 = re(F,/RC)*, 


1H. Scheffé and H. A. David have pointed out a serious error here. If there is an outlier 
at position (7, 7) then the assumption seems fair that s;2/o? is distributed as chi-square with 
(re — 1) degrees of freedom. If however there is no outlier then the mean-square of all rc 
residuals (MSD) is distributed as chi-square with rc degrees of freedom, and s;? is not chi- 
square distributed with (rc — 1) degrees of freedom. I do not now know how great are the 
effects of this error on the values given in Tables 1 and 2. The Editor has asked that the paper 
be published as given in the interests of having all papers and discussion in one issue. Re- 
vised tables will be published in this Journal as soon as possible. 
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we get the critical values given by the third member of that equation to compare 
with the statistic ¢, , the corresponding quantity obtained from data, and defined 
by equation 13. 


(13) t, rig max/8; ° 


Table 2 gives some approximate critical values of this ratio for an a per 
experiment of 0.2. 


TABLE 2 
Values of tig (a = 0.2) for R X C layouts. R and C, $ to 9. 


C 


| 


The statistic ¢, found from the data of Table 1 is 9/2.47 or 3.65, which is 
greater than the 3.48 given in Table 2. With a 20% risk of false alarm, then, 
the value at (3, 3) would be judged an outlier. 


3. Onze Brasep VALUE IN A 2” UNREPLICATED EXPERIMENT 


A further complication arises here in that the regression equation is to be 
derived from the data and is not specified beforehand as it is in the two-way 
layout. 

Routine analysis of the unrevised data may not suffice to give clues to the 
real effects. This will generally be the case if a defective value which has an 
expected value near one end of the range of responses actually appears near 
the. other end of the range. I see no alternative in such cases to making a series 
of trials, replacing in turn each of the 2” responses by the average of the others, 
and then revising the estimates of all effects accordingly. If one and only one 
of the resultant sets of effects permits a simple regression equation, containing 
mostly main effects and interactions among the factors having large effects, 
then a strong case can be made for naming as defective the value whose revision 
caused the simplification. If none of the substitutions proposed produces a 
notable simplification, two alternatives suggest themselves. There may be two 
bad values, which can be searched for using clues given by the one-at-a-time 
analysis just described. This failing, a half or full replication of the experiment 
may be necessary. But these grim contingencies are relatively rare. More common 
is the appearance of a single bad value, large enough to cause trouble, but 
not so large or so placed as to conceal all real effects. 

If a bias of 16 were to hit the treatment combination abed in a 2* experiment, 
the expected pattern of residuals, regressing only on the mean, is given by the 
first column of numbers in Table 3. If a different observation, not abcd, were 
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TABLE 3 
Effects on residuals of a bias of 16 at abcd in a 2 exp. 


Successive regression equations 


= ae i ee 
B c AB D ABCD 


| 
CSCOSCROOCOONNWND 
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_ 
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too large by 16, the residual 15 would of course shift to the corresponding row, 
all other rows showing —1. If, still assuming one value too large by 16 at abcd, 
we had found that effect A was large, we would have computed discrepances 
from a regression equation that contained an estimated mean and an estimated 
A-effect. The expected values of the resulting biased residuals are given in 
Table 3 in the column headed A. If the mean, and A and B were judged real, 
the discrepances would have the biases appearing in the following column. 
Similar remarks hold for the sequence of regression equations obtained by 
successively including the terms at the heads of the columns. 

The magnitude of the largest bias decreases with (16 — K), where K is the 
number of terms in the regression equation. Thus it might be thought that 
the inclusion of more terms will give a better fit to data that contains one bad 
value. However, the mean square for lack of fit does not continue to decrease, 
but settles down to a nearly constant value as indicated by the limiting case 
exhibited in Table 3, where the random errors are all zero. This persistent 
failure of the successive equations to fit better may well be a good indicator 
of the presence of a bad value. It would noi be an assurance that our mean 
square was converging on the true variance. 

In a 2” experiment, a bad value enters into every contrast, either increasing 
or decreasing its value by the amount of its bias. The average absolute value 
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of the contrasts will be increased, and so a half-normal plot (1) of the ordered 
contrasts will have too few contrasts near zero. A straight line through most 
of these contrasts does not go through the origin as it should but hits the abscissa 
axis at a positive value. This graph can be used, then, even before any dis- 
crepances are plotted, to warn us of a single biased value. It has the advantage 
that the detection does not depend on which effects are judged to be real. 

The pattern of biases depends on the location of the biased value and on 
the regression equation chosen to represent the data. The pattern can be quickly 
written down for any combination of interest with the aid of the standard 
transformation matrix of plus and minus ones that specifies how each observa- 
tion enters each effect-contrast. Such a matrix is given in Table M of the book, 
Design and Analysis of Industrial Experiments, edited by O. L. Davies. It is 
only necessary to mark the line corresponding to the suspected value and the 
columns corresponding to the desired regression equation. After reversing the 
signs of all entries in the columns in which the bias enters negatively, it is only 
required to add the plus and minus ones in each row to get the expected values 
of the propagated bias. Subtracting these values from the actual bias (a column 
of zeroes and one 16) gives the biases in the discrepances which are the c; listed 
in Table 3 and used in equations 14 and 15 below. 

Exactly analogous to the maverick contrast for two-way layouts, one can 
be written down for the 2” experiments. Thus 


2N 
(14) 6, _ DD ¢: d; — 2" Ge 


The variance of this contrast is 


7 (15) V(6,.) =o > = 2°(2" — ko’ 


t=1 


By reasoning parallel to that outlined above in equations 7 through 13, the 
critical ratio 


(16) the = (2"*/2"”)Fi 


can be derived. Table 4 gives some approximate values for this ratio for experi- 


ments with 16, 32 and 64 run factors and withf2,{4, 6 or 8-term regression 
equations. 


TaBLe 4 
Values of th, « (a = 0.2) for 2" experiments; K terms in regression 


n 4 5 6 


* See previous footnote. 





156 : CUTHBERT DANIEL 


REFERENCES 


[1] Danrst, C., Use of Half-Normal Plots in Interpreting Factorial Two-Level Experiments, 
Technometrics 1, No. 4. 

[2] Danret, C., Fractional Replication in Industrial Experimentation, Proceedings National 
Convention, American Society for Quality Control, 1957. 

[3] Dyxsrra, O., Partial Duplication, Technometrics 1, 63-75. 

[4] Scuurrt, H., The Analysis of Variance, Wiley and Sons, New York, 1959. 





TECHNOMETRICS May, 1960 


Discussion of the Papers of 


Messrs. Anscombe and Daniel 


Editor’s note: The two preceeding papers in this issue of Technometrics: ‘‘Rejection 
of Outliers” by F. J. Anscombe and “Locating Outliers in Factorial Experiments’ by 
C. Daniel, were presented at the Washington D. C. meetings of the American Statis- 
tical Association on December 29, 1959 at a session on the general topic of “Outliers.” 
The session was part of the program organized by the Section of the Physical and Engi- 
neering Sciences of the ASA. The papers were discussed by W. Kruskal, T. S. Ferguson 
an! J.W.Tukey. At the conclusion of the meetings the Editor proposed that the papers 
and their discussion be published in T'echnometrics and all participants agreed. Pro- 
fessor E. J. Gumbel, who was unable to speak at the meetings because of limitations 
of time, was also invited to present his comments. In a subsequent study, Mr. Daniel 
discovered an error in his paper that had been missed by the reviewers. The error is 
noted in Mr. Daniels paper in this issue of Technometrics. 


WiLuiAM KrusKAL 
University of Chicago 


The two papers of this session are obviously outliers among the population of 
papers appearing at statistical meetings. Further, these outliers are clearly in 
the positive direction. I have established this by use of a subjective rule, which 
considers how much papers stimulate me to think about the problems they 
treat and how much new ways of thinking about the problems are suggested. 

Both of the papers depart from the conventional treatment of outliers in 
that more complex situations than simple random normal samples are con- 
templated. Specifically, both papers discuss outliers primarily in the context 
of Model I analysis of variance designs with a high degree of symmetry. Both 
papers deal with rules for rejecting apparent outliers. The two papers differ 
in that Daniel’s is mainly concerned with rejection rate and power, while 
Anscombe’s is mainly concerned with. the mean square error of estimators 
following application of a rejection rule. I think that both papers help towards 
closing the ubiquitous gap between practice and theory. 

As to practice, I suggest that it is of great importance to preach the doctrine 
that apparent outliers should always be reported, even when one feels that 
their causes are known or when one rejects them for whatever good rule or 
reason. The immediate pressures of practical statistical analysis are almost 
uniformly in the direction of suppressing announcement of observations that 
do not fit the pattern; we must maintain a strong sea-wall against these pres- 
sures. The extra work and publication space needed to describe apparent out- 
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liers and the reasoning behind their treatment is usually negligibly small. 

Second, still about practice, one of the most important results of finding an 
outlier can be the correction of a flaw in the measurement process, or—even 
better—the creation of new insights into the phenomena under study. For 
such purposes, it may be that something like significance testing is more germane 
than evaluation of mean square error. 

In this general connection, I should like to call attention to a debate between 
H. C. Hamaker, on the one hand, and R. C. Anderson and A. W. Wortham, 
on the other. This appeared in the July, 1959, issue of Industrial Quality Control. 

I have two comments about theory. First, there may be cases in which the 
probability of a blunder depends on the value that would have been observed 
if the blunder were not present. For example, a scale reading may be more 
subject to gross misreading if it is towards the end of the scale than if it is 
near the middle. Can reasonably realistic models for this be stated and worked 
on? This question is closely related to another that I think deserves more 
attention: suppose that the probability that an observation is missing depends 
upon the value that would have been observed; to what extent are we justified 
in applying traditional missing-value manipulations? 

My final comment is more technical. In both papers, rejection rules are 
contemplated that sometimes lead to analyses with high symmetry except for 
one or two missing observations, missing now because of rejection by some 
rule. We know how to handle such situations with relative ease by computing 
the least-squares estimate of the missing observation’s expectation, putting the 
result in the sample in lieu of the missing observation, and then proceeding 
with estimation as if the full sample were there. These all-or-nothing rejection 
rules may have been proposed because of the ease of analysis following rejection 
of one or two apparent outliers. 

I suggest that alternative rules are worth consideration, in which the apparent 
outlier is not banished from the company of decent observations, but rather 
is given a lower weight. This can be done with hardly any more trouble than 
finding the so-called replacement for the missing value that one would do after 
rejection. 

Specifically, suppose that the sample is y, , --- , y, , with mean specified 
to lie on a given linear manifold in n-space. Suppose, to begin with, that 
Yr, *** » Yn-1 all have variance 1, but that y, has variance o”, a known number. 
(It would be enough if only the ratio between o” and the common variance 
of yi: , *** , Ya-1 Were known.) Let y be the estimate of the expected value of y, 
that one obtains following the usual missing-value approach for y, . Then one 
may state a linear combination of y and y, , (y + Ayn)/(1 + A), such that, 
if you replace y, by this linear combination and apply traditional analysis 
to the altered sample, you get the appropriate least-squares estimators with 
proper weighting for y, . The weight X is readily obtained in terms of the quan- 
tities needed to find y. 

This may be derived algebraically by the use of M. A. Woodbury’s result 
on inverting special matrices, or it may be obtained geometrically. Extensions 
to more than one separately weighted y,; are undoubtedly possible in various 
ways, and I hope to return to this topic in more detail elsewhere. 
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DISCUSSION OF ANSCOMBE AND DANIEL PAPERS 


Tuomas S. FERGUSON 
U. C. L. A. and Princeton University 


I wish to thank and compliment the speakers for producing first rate work 
on the theory of rejection of outliers. Both papers treat a topic which has not 
even been attempted before in the outlier literature, namely, the problem of 
rejection of outliers in designed experiments, where the means of the observa- 
tions may be complicated functions of unknown parameters so that even dis- 
covery of a bad value may be difficult. The methods of the two papers differ. 
Whereas Daniel makes use of “‘patterns” which appear in the residuals, Anscombe 
goes straight to the extreme residual for his rejection statistic. The fact that 
the correlation test that Daniel proposes turns out to be equivalent to the 
studentized extreme residual is both surprising and gratifying in that it is a 
studentized justification of the rejection rule of Anscombe, at least in the two 
designs considered. One suspects that this could be extended at least to all 
designs considered by Anscombe. Another justification can probably be made 
by extending the result of Paulson (“An optimum solution to the k-sample 
slippage problem for the normal distribution,’ Annals of Math. Stat. 1952,) 
to two-sided rules and for designed experiments. 

Daniel’s treatment of the two-way layout is a nice neat package. I want 
only to say that I don’t think one can avoid using a multiple comparison method 
here. The extreme residual is not a simple contrast, but, however, it is a function 
of all the simple contrasts. 

Finding outliers in factorial experiments is inherently more difficult because 
one doesn’t know which effects are going to have large mean values. Further- 
more, one maverick may mask the real effects so that we may not be able to 
judge which effects to take into the regression equation. The use of the half- 
normal plot here is ingenious. However, one needs a more satisfactory way 
of finding which of the original observations is the bad one. The original half- 
normal plot may tell you that there is a spurious observation but one would 
hate in a 2° experiment to have to make 32 more plots, one each for treating 
each original observation in turn as a missing value, in order to determine 
which observation is the bad one. 

A related point is of interest here: that of chosing the plotting points for the 
half-normal plot. Suppose we have an ordered sample of size n from a half- 
normal distribution z, < x. < --- < 2, . For the half-normal plot we plot 
the points (x; , y;) in the plane, where the y; are chosen beforehand so that 
the points will in some sense lie on a line. It is usual to take y; as that number 
such that F(y;) = 7/(n + 1) (or F(y;) = (¢ — 4)/n) where F is the distribution 
function of the half-normal distribution. I would prefer to use y; either as the 
expected value of x; , or as the ith entry of the eigenvector corresponding to 
the largest eigenvalue of the second moment matrix of the x; . The first has 
the property that x-intercept of the least squares regression line of x on y is 
an unbiased estimate of zero. The second has in addition to this the property 
that the expected value of the residual sum of squares of the distances to this 
least squares regression line is a minimum. However, tables of the expected 
values and of the second moments of the order statistics for the half-normal 
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distribution do not exist. I would like to appeal to the statistical public to 
compute such tables. 

I feel that Anscombe’s paper will be the cornerstone on which tatiies genera- 
tions of outlier rejectors will base their work. The basic idea, which appears 
is all statistics, is that one can better design the experiment if one knows what 
actions one will take later. This is a point which decision theory emphasizes. 
However the proferred solution to this problem took an extra bit of ingenuity; 
it is not decision theory but insurance. I would like to make only one remark, 
that Rule 0, studentized, is exactly Thompson’s rule, (Annals of Math. Stat. 1935) 
and does not suffer from the defect mentioned in the paper for the nonstudentized. 
I think it is a serious contender for attention in the extension of these ideas 
to the problem when the variance is unknown. 


Joun W. TuKey* 
Princeton University and Bell Telephone Laboratories, Murray Hill 


I would like to begin by saying that I agree with all the good things that the 
earlier discussants have said about the papers. 

This morning’s meeting illustrates how important it is for us to return to 
the “Elementary Problems of Statistics”. One of the previous discussants sug- 
gested that the two papers would be foundation stones for the future work 
on the rejections of outliers. Notice that these papers were not thought of as 
finishing a process but rather as a beginning, as initial steps in further growth 
of work on outliers. Surely the rejection of outliers is one of the really elementary 
problems of statistics. But, as we are coming to recognize, it does not represent 
a problem which is anywhere nearly completely treated. 

In connection with the discussion of whether a rejection procedure should 
be a test of significance or not, I think that-it is important to emphasize the 
way in which Frank Anscombe’s paper avoided the semantics of “rejection’’. 
He was concerned with doing the best that he could with the data before him, 
with combining it in the best way that seemed reasonably feasible ... not 
with whether an individual observation should be dropped or should not be 
dropped, not with whether an individual observation was in some specific sense 
unusual or not ... only with the question of what was the best thing to do. 
By asking how well one does in estimation, ene avoids any implication of making 
up one’s mind whether or not a particular observation is a “bad value”, one 
avoids any need for “significance testing’. 

Bill Kruskal has pointed out the connection between such yes-or-no procedures 
as rejection and the much more flexible and quantitative procedures of weighting. 
(At this point it may be appropriate to notice one further related reference, 
a paper by H. R. Hulme and L. S. P. Symms called ‘“The Law of Error and 
the Combination of Observations’, Monthly Notices of the Royal Astronomical 
Society, June 1939, Vol. 99, pp. 642-649. This paper discussed the use of a 
weighting procedure in reducing some 9,000 observations made with a floating 
zenith telescope at Greenwich over a ten-year period. It is, of course, related 


* Prepared in part in connection with research at Princeton University under Contract No. 
DA 36-034-ORD-2297 sponsored by the Office of Ordnance Research. 
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to the Jeffreys techniques referred to in Frank Anscombe’s paper.) There are 
three qualitative operations which can be used to deal with outliers and wild 
values. It may well be that the best practical solutions will turn out to involve 
combinations of these. These are (i) rejection, (ii) truncation, and (iii) Winsoriza- 
tion. We have heard a lot today about rejection, where one chooses which 
extreme values to delete on the basis of the specific and quantitative relation 
of these values to the other values observed. In symmetrical truncation, by 
contrast, one deletes a certain percentage of high values and an equal percentage 
of low values, without regard to just how the values compare with one another. 

Finally there is Winsorization. Charles Winsor put forward a principle of 
quite general application, namely: While the numerical value of an apparently 
wild observation is untrustworthy, the direction of its deviation (e.g. high 
or low) is worthy of some attention. He applied this principle to outliers by 
taking the largest deviations or largest residuals and decreasing their magnitude, 
while retaining their sign, until they are equal to the next largest ones, thus 
making a qualitatively reasonable adjustment. In fact, it is possible to give 
quite quantitative reasons why this sort of adjustment is not only a convenient 
approach but an effective one. 

The natural way to attack the class of problems of which defense against 
wild values js a part, is not to begin in small samples with rejection, but to begin 
in large samples with Winsorization and truncation. However one may, and 
I believe should, look forward to an eventual tendency to use a combination 
of two or more of these procedures in all situations. The tentative suggestion 
in Cuthbert Daniel’s paper that, for example, in a 2° factorial experiment one 
might plan to reject one observation all the time and a second one some of the 
time, is precisely a proposal to combine some truncation with some rejection. 

(I would like to say that I agree with him both in putting forward this idea 
and very particularly, in putting it forward tentatively. It seems to me this 
is the sort of procedure that must be considered, but it is also clear that it is 
the sort of procedure which, under some circumstances, might give rise to 
difficulty. We will have to think about this sort of procedure carefully, and 
try it out in practice, before we know when and where it can be routinely 
adopted.) 

I would like now to say a few things about the case where o is known. The 
results in Anscombe’s Table 3 can be put into somewhat more compact form 
if one notices that, for fixed premium and to an adequate approximation, C is 
proportional to (v/n)°’*, a is proportional to (x/n)~***, and b is proportional 
to (v/n)***. If we take advantage of such facts we can greatly compress the 
tables upon which we are to base our judgments. 

Frank Anscombe emphasized the size of the premium that had to be paid, 
and he also emphasized that the amount of protection provided had also to 
be considered. I suspect that the protection here is perhaps best thought of 
in terms analogous to the collison insurance of automobiles, in terms of “what 
is the amount deductible; what is the largest net loss which one may have to 
suffer?” It would clearly be of interest to know just how much of an increase in 
variance one is going to suffer in the worst situation. 

It is important to realize that there are many different sorts of difficulty. 
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The particular sort of difficulty that has been emphasized in Anscombe’s paper 
is exactly that which arises when a normal distribution is supplemented by 
two isolated masses of probability, one at +bo and one at —bc, and these 
isolated masses are just big enough so that one observation drawn from them 
will occur in at least a moderate fraction of samples of the chosen size, but 
two observations per sample will occur quite infrequently. Clearly one must 
think about different values of b. It even may be that one can think about a 
probability distribution for these values of b. But it is doubtful that any one 
probability distribution would suffice. Throughout the discussion of outliers 
and wild values we have got to face a variety, all tentative, of bad possibilities. 

As some of you in the audience know, I tend to consider cases where the 
badness is represented by the mixture of two normal distributions, with the 
same mean but with different variances (where we have a contaminated distri- 
bution), much more than the case where we have a normal distribution plus 
some specific values with finite probability. Here, too, one must consider dif- 
ferent possibilities, exemplified by different degrees of contamination. Let us 
look at the case where we have a large sample and are worried about this other 
sort of difficulty. If we consider contaminated distributions in which the mixture 
includes a certain fraction y of a wider normal distribution with three times 
the scale of the narrower norma] distribution (which appears at a fraction 1 — y), 
then we have the results shown in the figure (taken from ‘‘A Survey of Sampling 
from Contaminated Distributions”, Chapter 39 of Contributions to Probability 
and Statistics; A volume dedicated to Harold Hotelling, Stanford University 
Press, 1960) for various sorts of truncated means. 

The curve labeled 0%-truncated is of course the curve for the ordinary mean, 
where nothing has been deleted from the observations. At the bottom we have 
the curve where we have removed the 50% highest observations and the 50% 
lowest observations, leaving only a single observation in the middle, so that 
the mean of the retained values turns out to be mean of a single central value, 
that is the median. The other curves show the result of cutting one, three, 
or six percent, respectively, of the observations off the upper end of the sample 
and an equal amount off the lower end of the sample. These results are expressed 
in terms of asymptotic variance and large sample theory; they are surely not 
applicable down to samples as small as say 10, though they may be fairly effective 
for samples of the order of 100. Surely, they do not work down to the level 
of samples of three or four, which Frank Anscombe pointed out are of frequent 
occurrence. Here the premium (in Frank Anscombe’s terms) is the amount 
by which the curve lies below 100% for y = .00. For as much as six percent 
truncation on each tail, the premium is only about 3%. For a 3% premium 
we can protect ourselves against a possible loss of about 30%. This example 
makes clear the possibility of using truncation in large samples. 

Next a few words about the case where o is unknown, the case where we are 
concerned either with regresson situations regarded as such, or with patterned 
experiment situations which can be treated from a regression viewpoint. While 
Cuthbert Daniel was clearly aware that there was no need to calculate the 
standard deviation of the cleaned-up observations in order to get something 
to compare with the maximum residual, I should like to emphasize that this 





DISCUSSION OF ANSCOMBE AND DANIEL PAPERS 


6% TRUNCATED 


3% TRUNCATED 


1% TRUNCATED 


"THE MEAN" 
=0% TRUNCATED 


> 
oO 
2 
¥ 
S 
w 
w 
uJ 
Oo 
Ee 
Oo 
- 
a 
= 
> 
” 
2 


"THE MEDIAN"= 50% TRUNCATED 


00 Ol 02 .05 
FRACTION OF CONTAMINATION, Y 


Asymptotic efficiency, for location, of truncated means. 


(Figure reprinted with permission from Contributions to Probability and Statistics, Olkin et al., 
Copyright 1960, Stanford University Press.) 


is the case. I should also like to emphasize that one can always compare the 
maximum resid’1al with a suitable multiple of either the square root of the sum 
of squares of residuals, or with a different suitable multiple of the square root 
of the mean square for residuals. The actual numerical value of either suitable 
multiple will depend, in a two-way case, for example, upon the number of 
rows and columns used. (Such dependence is going to occur for the critical 
value of any ratio one chooses in this situation.) One can always compare the 
maximum residual with something calculated from all the residuals together. 

I cannot refrain from joining the authors and earlier discussants in emphasizing 
again that all the most natural ways of dealing with the residuals tend to be 
mathematically equivalent; that using the correlation coefficient, or the relatively 
biggest residual, leads to the same basis for making choices whether or not to 
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reject. And I also cannot refrain from underlining the fact that, in the long 
run, we are going to need to know somewhat more about the joint sampling 
behavior, in a wide variety of situations, of residuals. 

Comments, both by: speakers and discussants, make it, I think, also necessary 
to say a word or two about some multiple comparisons questions. It has been 
suggested that it is necessary to use Scheffé’s approach. It seems to me that 
it is only necessary to make some sort of allowance for the fact that, when we 
look at the largest residual, we are looking at the largest of a certain number 
of possibilities (which have a certain interrelation), that we are doing some 
selection. If one were to use the distribution of this largest residual itself to 
generate a, multiple comparisons procedure, one would get a somewhat more 
stringent ‘‘5% level” than if one used an F-statistic-based multiple comparison 
procedure at a ‘5% level”. And the latter, again, would be more stringent 
than a procedure in which, when there were R rows and C columns, we would 
say that we have RC residuals, so we will take a critical value of |¢| which cor- 
responds not to 5% but to (5/RC)%. There will be some differences between 
these three approaches, but the differences will be minor and unimportant. 
The important thing is that all allow for the problem of selection. Neither the 
detailed approach nor the precision of its application will be important. 

There was a suggestion by the discussants that Cuthbert Daniel’s paper 
was concerned with model 1 anova (model 1 analysis of variance). I have to 
point out that this is not really true, especially because one of the things with 
which he was explicitly concerned in his talk was the fact that the particular 
sort of regression, the particular sort of separation between accepted effects 
and residuals, which would be appropriate for a given set of data, would depend 
upon what the data looked like after one (or more) of the values had been 
deleted. To have flexibility in the way in which the data are going to be inter- 
preted, flexibility which is directly associated with the data values themselves 
is not to have model 1 anova. It is to have something else, something which 
is harder to analyze from the point of view of outliers, something which is much 
more reasonable from the point of view of actual data analysis. One of these 
days, perhaps, it will have a specific name, but for the present it should suffice 
to point out that it is not model 1. 

There are a few other specific points that might perhaps be mentioned. One 
is that, if (i) the experiments are performed one after another and (ii) blocking 
is not too serious a problem, then, rather than planning to reject the experiment 
with the largest residual and analyze the others, one might plan to keep a con- 
tinuing check on how the results are coming out, so as to know in advance of 
the completion of all the experimental runs (i.e. experimental trials, or plots), 
which run or runs seem likely to give the largest residuals so as to be able to 
set aside these runs providing large residuals, and to repeat each of the cor- 
responding experimental conditions twice as many times as before. This would 
give good basis for deciding whether such wild shots occur randomly, or whether 
they are really associated with certain combinations of treatments. This is 
basically one of the most important problems in connection with outliers, one 
which has already been stressed by Bill Kruskal. It is often very hard to deter- 
mine whether one wants to set the apparent wild shot aside and produce a 
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good analysis of what is left, or whether, since such understanding will be more 
valuable than any possible increase in the accuracy of the pre-chosen estimate 
or estimates, one wants to cling to the wild shot above all else as something 
whose occurrence must be understood. 

If, on the other hand, the experiments are performed in parallel time-wise, 
so that duplication of the post facto unusual ones is out of the question and 
only deletion is possible, this problem: ‘Which is the baby, and which is the 
bath-water?” becomes, of course, very much more stringent. In either circum- 
stances, whenever the wild shots (those things which appear from a factorial 
point to be “‘wild shots’) are such as to indicate that some other way of looking 
at the data analysis would be better, one must be prepared to give up the 
factorial approach. 

At one time Cuthbert Daniel objected to the use of a log-log type of half- 
normal plot because it made the smallest effects too noticeable. Now he is 
concerned with values moderately far down in the size ordering of effects. 
Now that such values are used in judging whether or not there appears to be 
a wild shot present, there will perhaps be more interest in the log-log form of 
the half-normal plot. 

Finally, in closing, I would like to ask a question of interest, but probably 
not immediately of serious import. When are we going to get a combination 
of the problems we discussed at the Pittsburgh meeting and the problems we 
discussed here? When are we going to have the first discussion of the rejection 
- of outliers in supersaturated patterns where the value of v/n is already negative? 


E. J. GuMBEL 


Columbia University 


The acceptance or rejection of outliers is an important decision problem for 
all observations. Thus, this question arises for any statistical distribution used 
in the analysis of numerical data. 

The outliers are values which seem either too large or too small as compared 
to the rest of the observations. Thus, they are extremes. The difficulty of the 
decision arises from the fact that the absolute size of the extremes increases in 
general with the sample size. A value which may be rejected for a certain number 
of observations has to be accepted for a larger sample since the gap between 
the outlier and the rest will be filled by new observations, one or several of 
which may even exceed the old outlier. Therefore, any criterion for such a 
decision must be based on the sample size. Anscombe proposes, for a normal 
population with known o, to reject the absolute largest among n observations 7, if 

lan a pl > Co, 
where yu is estimated as the mean of the retained observations. If o is unknown 
it is to be replaced by s. 

Anscombe’s notation does not bring out the fundamental fact that C increases 
with the sample size n. Therefore it should be written C,, . A table of C, as a 
function of n for chosen values of a would be desirable. 

Anscombe claims that for n = 3 and 4 observations assumed to be normal with 
a known standard deviation “this method is probably the best that can be 
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suggested.”’ The validity of this statement depends of course on what we define 
as being the best. Logical consistency, numerical simplicity of procedure and other 
criteria may lead to different definitions for such an optimum. 

All rejection criteria are based on the extreme values as linked to the sample 
size. The simplest procedure for the normal distribution seems to me to start 
from the distribution of the extreme normal values for n observations and to 
choose an appropriate probability a for the rejection. The choice of this a is 
arbitrary and the conventional value of .05 may be appropriate. This method 
can be generalized for any known initial distribution. If it is not symmetrical the 
two distributions of the largest and of the smallest values as functions of the 
sample size have to be considered separately. Then Anscombe’s rules could be 
applied in a modified way to these distributions of the extremes. 

The rejection of outliers on a purely statistical basis is and remains a dangerous 
procedure. Its very existence may be a proof that the underlying population is, 
in reality, not what it was assumed to be. The alleged homogeneity may not be 
true and the outlier may be the most important observation because it throws 
light on the way the observations were made. 

In the history of this problem, the important contribution of P. R. Rider, 
“Criteria for rejection of observations’, Washington University Studies No. 
8,1933 and the literature given therein should be mentioned. 


AUTHOR’S REPLY TO DISCUSSION 


F. J. ANSCOMBE 


This has been a valuable discussion, because a number of suggestions have 


been made that wil! be really helpful in further research, and also because it 
serves to set the two papers in a just perspective. There are so many matters 
that can be discussed regarding outliers, so many different situations in which 
a problem of outliers arises, that one must be careful to note just what topic 
or situation is under consideration at any point. The topics that Cuthbert 
Daniel and I have discussed are narrow and special, relative to the whole field 
of outlier debate—though they may have seemed broad enough while we worked 
on them. 

In comment on Cuthbert Daniel’s paper, I should like to express the hope 
that his idea of detecting outliers by examining half-normal plots will be pursued 
further. This appears to be an essentially different device from those based 
on the largest residual. 

I hope that further work concerning outliers and residuals will be presented 
at the Fourth Berkeley Symposium. 

In revising my paper for publication, I have made various small amendments, 
prompted by the above public discussion, by private discussion with Churchill 
Eisenhart, and by detailed and searching comments from Leonard J. Savage 
and William Kruskal. To all I offer my grateful thanks. Irwin Guttman and 
I regret that the program of computations undertaken last year was not quite 
completed. We hope to add a supplement to Table 4 in a later issue of this 
journal. 
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Tests for the Validity of the Assumption that 
the Underlying Distribution of Life is 
Exponential: Part II 


BENJAMIN EPSTEIN 


Stanford University and Wayne State University 


Part I of this paper was published in the previous issue of Technometrics, (Vol. 2, 
No. 1, February 1960). In Part I Dr. Epstein describes several graphical and 
analytical procedures for testing the assumption that the underlying distribution 


of life is exponential. This portion of the paper contains numerical examples illustra- 
tive of these procedures. 


INTRODUCTION 


In Part I of this paper the following tests, useful to determine the validity 


of the assumption that the underlying distribution of life is exponential, are 
described: 


. A Graphical Procedure 

. The x’ Test of Goodness of Fit 

. A Criterion Based on the Conditional Distribution of Total Lives 

. A Test for Abnormally Early Failures 

. A Test for an Abnormally Long First Failure 

. A Test whether the Mean Life (or Rate of Failure) in the First Half 
of a Life Test Differs Significantly from the Mean Life (or Rate of 
Failure) in the Second Half of a Life Test 

. A Test for Whether or Not the Mean Life Fluctuates During the Life Test 

. A Special Procedure for Procedure (7) When r = 1 

. A Test Based on the Maximum F Distribution 

Tests for Abnormally Long Periods in which there is No Failure 
. A Graphical Procedure Based on the Kolmogorov-Smirnov Test 
. A Test Based on Conditional Rate of Failure 


The description of these tests is followed by general comments, a short summary, 
an appendix in which some of the more mathematical arguments are found, 


and finally by a bibliography. Part II of this paper provides numerical examples 
descriptive of the procedures given in Part I. 


EXAMPLES 


(1) 49 items were placed on life test and the test was run until all items failed. 
The observed failure times 7, , --- , T49 Were: ~ 


167 
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72.4 102.8 151.6 

73.6 108.5 152.6 

76.8 128.7 164.2 

‘ : 83.8 133 .6 166.8 

: . 95.1 144.1 178.6 
25.1 ; 97.9 147.6 185.2 
35.8 ; 99.6 150.6 187.1 


Test the hypothesis that the underlying distribution of life is exponential, 
using procedure 1. 


Solution: According to procedure 1, the graph of 7; against 


o n+1 ) 
=. toe (RL 
should be approximately a straight line passing through the origin, if the under- 
lying distribution is exponential. 

In our example we plot 7; against (n + 1)/(n + 1 — 7) on semi-log paper 
(see Figure 1). 


NNO Oe 
eOwMmOoo nN WY 


— 


(2) 49 items were placed on life test and the test was run until all items failed. 
The observed failure times 7, , «++ , T49 were: 


162 287 539 884 1178 1666 2181 
172 301 629 886 1205 1676 2193 
199 302 634 918 1377 1742 2445 
200 389 643 1008 1406 1818 2741 
218 393 682 1101 1501 1896 3191 
220 401 706 1129 1576 2002 3577 
271 508 777 1146 1641 2021 3694 


One suspects that the waiting time for the first failure was abnormally long, 
and that the underlying p.df. is 1/0 e*"*”", with A > 0. Give a graphical 
procedure for testing whether A > 0. 

Solution: Again, as in the previous example, we plot 7; against 


1 
v= oe (FF) 


As before, the data are plotted on semi-log paper (Figure 2). If the underlying 
distribution is one-parameter exponential, we should get a straight line passing 
through the origin (0, 0); if the underlying distribution is two-parameter expo- 
nential, we should get a straight line passing through (A, 0). In our example 
it is fairly evident that we would reject the hypothesis that the distribution 
is one-parameter exponential and would accept the hypothesis that it is two- 
parameter exponential (see Figure 2). 

Remark. One may also find an approximate estimate of A by looking at 
the graph. In our particular example, A is approximately equal to 150. 

In example 8, we give a numerical procedure for testing whether A > 0, 
based on procedure 5. 
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150 180 210 
OBSERVED FAILURE TIME 


Figure 1—Graph of y; = log (n + 1)/(n + 1 — 7) against 7; . 
Data from Example 1 (n = 49). 


OBSERVED FAILURE TIME 


Figure 2—Graph of y; = log (n + 1)/(n + 1 — #) against 1; . 
Data from Example 2 (n = 49). 
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(3) One observes a stochastic process until 7,; , the time when the 13th event 
occurs. The 13 observed times are r, = 5.6, 7, = 12.8, r3 = 17.6, 7, = 27.6, 
Ts = 44.2, 7, = 45.5, t7 = 68.8, rs = 89.2, tr, = 90.7, t19 = 100.2, 71, = 106.5, 
tig = 115.5, 715 = 123.8. 

Are these observations consistent with the hypothesis that the underlying 
process is Poisson with constant rate parameter? Use procedure 3. 

Solution: If the Poisson assumption holds, then the 12 times 7, , 72, -+* , 712 
should be tested to see whether they are uniformly distributed in (0, 7,3). The 
theoretical expectation of >>}2, 7; is (12/2)r1, = 6(123.8) = 742.8. The theo- 
retical standard deviation is (12/12) (123.8) = 123.8. Under the hypothesis 
that the underlying process is Poisson with constant rate parameter, >.}?, 7, 
is approximately normally distributed with mean 742.8 and standard deviation 
123.8. Hence a 95% acceptance interval for the hypothesis being tested is 
742.8 + (1.96)(123.8) = 742.8 + 242.6 = (500.2, 985.4). The observed sum is 

21 7; = 724.2. Since this falls within the acceptance interval, we accept 
the hypothesis that the underlying process is Poisson with constant rate param- 
eter at the .05 level of significance. 


(4) Twenty items are placed on test. The test is discontinued after 11 failures 
occur. Items are not replaced as they fail. The first 11 failure times.are 7, = 7, 
T, = 12, 73 = 15, tr, = 24, 75 = 25, 7. = 48, 77 = 53, 7, = 56, To = 72, Tio = 95, 
711 = 100. Are these data consistent with the hypothesis that the underlying 
distribution of life is exponential? Use procedure 3. 

Solution: Let us compute V(r;) for 7 = 1, 2, --- , 11. Since the number of 
items originally on test is nm = 20, V(7,) = 20(7) = 140. V(r.) = V(r) + 
V(r.z — 71) = 140 + 19(5) = 140 + 95 = 235. Proceeding in this way we get: 
V(r1) = 140, V(r.) = 235, V(r3) = 289, V(r.) = 442, V(rs) = 458, V(r.) = 803, 
V(rz) = 873, V(rs) = 912, V(r9) = 1104, V(r10) = 1857, V(71,) = 1407. Under 
the assumption of exponentiality, the 10 total lives V(r;), 7 = 1, 2, --- , 10, 
should be uniformly distributed in (0, V(7:,)). To test this we calculate the 
theoretical mean and theoretical standard deviation for >.:2, V(r;). These are 
(10/2)V(r1,) = (10/2)(1407) = 7035 and +/10/12 (1407) = (.913)(1407) = 
1285, respectively. Under the exponential assumption, >»:2, V(z,) is approxi- 
mately normally distributed with mean 7035 and standard deviation 1285. 
Hence a 95% acceptance interval for the hypothesis that the underlying distribu- 
tion is exponential is given by 7035 + 1.96(1285) = 7035 + 2518 = (4517,9553). 
But the observed sum is }>32, V(r;) = 6613. This is well within the acceptance 
interval and hence the hypothesis that the underlying distribution is exponential 
is accepted at the .05 significance level. 


(5) In a non-replacemeni situation, a life test is terminated at a preassigned 
total life of 1000 hours. The total lives T(r,), «++ , T'(745) are as follows: 


30 173 266 346 468 611 662 763 915 
63 178 281 392 486 616 692 773 
99 207 290 433 562 631 693 776 

102 244 310 436 563 651 693 807 

172 255 343 455 584 660 700 844 
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where 7(r;) = nr, and T(r,) = Dorit: + (n — r+ 1)t,,7 = 2, «++ , 45. 

Are these observations consistent with the hypothesis that the underlying 
process is Poisson with constant rate parameter? Use procedure 3. 

Solution: Using procedure 3, it follows that if the Poisson assumption holds, 
then the total lives T(r,), --- , T(745) are uniformly distributed in (0,1000). 
Since the sample is fairly large, the x” goodness of fit test may be applied to 
see whether there are abnormal deviations from the uniform distribution. 
Dividing the interval (0,1000) into 5 equal intervals, we have* 


Expected Number (e;) 


Observed Number (0;) | 


5 2 
2. vy lai—e) _ 26 _ 
x= dX +: “~— 2.89. 
Since we have divided the interval (0,1000) into 5 subintervals or classes, 
*_1 [(0; — e:)"]/e; is distributed as x7(4). But Pr(x?(4) > 2.89) lies between 
.5 and .7. Hence we accept the hypothesis that the underlying process is Poisson 
with constant rate parameter. 


(6) In a life test we obtain the first 10 failure times. We suspect that the first 
2 failures may have occurred unusually soon and should be rejected. The total 
life associated with the first two failures is T(r.) = 24 and T(rio0), the total 
life associated with all ten failures, is 600. Test whether 7, and 7, occur ab- 
normally soon. Carry out the test on the .05 significance level, using procedure 4. 

Solution: Let us substitute in (11), using a = .05. The lower 5% point of 
F(4,16) is F** = 1/5.84, since Pr(F(4,16) < 1/5.84) = .05. Hence we will 
reject if 


2(600 — 24) 
T (12) < (5.84)8 ’ 
ie, if 


576 
T (12) < 23.36 — 244. 


Since T(rz) = 24, we reject the null hypothesis on the .05 significance level 
and instead have evidence that the first two failures occur abnormally soon. 


* Note that we have used classes of equal probability rather than equal length. This removes 
the ambiguity of chi-square that was discussed in our description of procedure 2. For a dis- 
cussion of this point, see Gumbel [20] and [21, p. 28]. Even when using classes of equal proba- 
bility, one ambiguity remains. How many classes should be chosen? A reasonable rule is to 
divide the interval in such a way that the common value of the expected number e; > 5. 
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(7) Twenty items are placed on life test. The life test is discontinued after 
11 items have failed. Items that fail are not replaced. The first 11 failure times 
are 7, = 107, 7, = 112, 73; = 115, = 124, 75 = 125, tr, = 148, 7, = 158, 
ts = 156, tr, = 172, tro = 195, and 7,, = 200. The first failure time appears 
to be abnormally long and one wishes to test whether the underlying distri- 
bution is (1/0)e"*"“”", with A = 0 against the alternative that A > 0. Carry 
out the test at the .01 significance level, using procedure 5. 

Solution: It can be computed readily that 7(7,) = 2140 and T(r,,) = 3407. 
Hence T(71, — 7:1) = 3407 — 2140 = 1267. Substituting in (13) with a = .01, 
we see that we will reject the hypothesis that A = 0, if T(7,) > (5.85) (1267)/10 = 
741 [5.85 is the upper 1% point of the F(2,20) distribution]. 

But 7'(r,) was observed to equal 2140. Hence we reject the hypothesis that 
A = 0 at the .01 significance level. We have strong grounds for believing that 
the minimum life A > 0. 


(8) Test whether A > 0, ie., whether the first failure time is abnormally 
long in example 2. Carry out the test on the .01 significance level using pro- 
cedure 5. 

Solution: In this problem TJ(7,) = 7938 and T(rts) = 58432. Hence 
T(t19 — 7:1) = 58482 — 7938 = 50,494. We will reject the hypothesis that 
A = 0 on the .01 significance level if 


a 

T(r.) > Gt ee =) 

where G*,, = 4.83 is the upper 1% point of the F(2, 96) distribution. Since 
the right-hand side is 5050 and 7(7,) = 7938, we reject the hypothesis that 
A = 0 at the 1% significance level. There is good reason to suspect that it 
took abnormally long for the first failure to occur. 


(9) A life test is carried out and one is given the first twenty failure times 
T ST. S +++ < too. The total life T(r,9) associated with the time interval 
[0, 710] is 10,000 hours and the total life T'(r2. — 710) associated with the time 
interval [7,9 , 720] is 25,000 hours. Let us define 0, as the mean life in [0, t10] 
and 6, as the mean life in [719 , 720]. (a) Test the hypothesis 0, = @, against the 
two-sided alternative that 6, ~ 0, , at the .05 significance level; (b) test the 
hypothesis that 6, = 6, against the one-sided alternative that 0, > 06, , at the 
.05 significance level. Use procedure 6. 

Solution: Using procedure 6, let us note that 27'(r10)/0, and 27729 — rT10)/02 
are mutually independent and are each distributed as x’(20). If, in particular, 
the null hypothesis is satisfied, 0, = @. = @ and so T(T20 — T10)/T(r10) is distri- 
buted as F(20, 20). 

But Pr(1/2.46 < F(20, 20) < 2.46) = .95. Hence if we wish to have a two- 
sided test with type 1 error = .05, we should accept the hypothesis that 6, = 4 if 


es 


T(T20 aan T10) 
2.46 oc 


< 
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and reject otherwise. The observed value of 


T(T20 = T10) ad 25 ,000 — 
T(t10) i 10,000 zs 


Hence we reject the hypothesis that #, = 0, on the .05 level. 
In the one-sided case, we note that Pr[F(20, 20) > 2.12] = .05. Hence we 
will reject the hypothesis that 6, = 6, (and accept 6, > @,) if 
T(t20 — 7 10) 
T(r.) > 2.12. 
Thus we reject the hypothesis that @, = 6, on the .05 level, since we observed 
T(t20 — T10)/(T10) = 2.5. 

Remark: In actual practice there may well be a prior grounds for using either 
the one- or two-sided test. Suppose we know, e.g., that the quality (as measured 
by length of life) of the last 10 items tested is at least as good as the quality 
of the first 10 items. If that is the case, a one-sided test should be used. 


2.50. 


(10) One observes a life test until 7;, , the time when the 51st failure occurs. 
The observed total lives are: 


150 1220 1961 2476 = =3270) =. 3575 3790 
314 1273 =©2165 + = 2789 3297 3577 3799 
493 1329 2180 2838 $3348 3577 3810 
517 1401 2275 2853 3357 = 3579 

858 1449 2282 §=©3006 «=6.: 3402—t—«é88'6211 

865 1549 2434 3191 3461 3656 

889 1716 2439 3228 3501 3694 
1033 1728 82441 3238 3509 3760 


It is suspected that as time goes on, there is a deterioration in quality so that 
the mean life 6 decreases as a function of time. Hence we are interested in testing 
the hypothesis that the observations are coming from an exponential distri- 
bution with constant parameter @ against the alternative that the parameter @ 
decreases during the life test. Use procedure 6 at the .05 significance level. 

Solution: In using procedure 6 we shall be testing a somewhat simpler hypo- 
thesis. Let us define 6, as the mean life in [0, 7.5] and 6. as the mean life in 
[72s , Ts1]. Then we should like to test the null hypothesis: 6, = 6, against the 
alternative: 6. < 6, . Using procedure 6, 27 (725)/0, and 27'(75,; — 725)/02 are 
mutually independent and distributed as x7(50) and x*(52), respectively. If 
the null hypothesis is satisfied, 6, = 0, = 6 and the statistic [257 (75, — t25)]/ 
26T (725) is distributed as F(52, 50). The null hypothesis will be rejected for 
sufficiently small values of F. Using the given data, we have 


257'(t51 — 725) _ 25(3810 — 2476) 
267 (725) 26(2476) 


From the F table we have that Pr[F(52, 50) < 1/1.60 = .625] = .05. But 
.018 < .625, hence we reject the hypothesis that 6, = 6, at the .05 significance 


= .518. 
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level and assert that 6, , the mean life in the second half of the experiment, 
is less than 6, , the mean life in the first half of the experiment. 

Remark: A graphical way of seeing whether 6 changes in time is to plot 
T (ri: — 7;), the elapsed life between successive failures, against 7. An analytical 
procedure is to compute a rank correlation coefficient for the {T(r;., — 7;)}. 


(11) We are given the results of 5 life tests in each of which 10 failures were 
observed. The total lives in each of the life tests are 12.74, 10.02, 7.25, 9.18, 
and 12.99 (thousand hours), respectively. Use procedures 7 or 9 and the formulae 
in Appendix 2 to test whether @ is constant over all the tests. 

Solution: Let us carry out the solution using Appendix 2. In this problem r, 
the number of failures in each group, is 10 and k, the number of groups, is 5. 
It is readily computed that 


log L, = {log 12.74 + log 10.02 + log 7.25 + log 9.18 + log 12.99] 


\ (52:18) 


= 2.3223 — 2.3453 = — .0236 
“. L, = .9772. 


Reference to the L, tables for n = 2r + 1 = 21 and k = 5, indicates that 
the hypothesis that 6 is constant should be rejected if L, < .905, if we wish 
to maintain a Type I error = .05. The value actually obtained is .9772 and 
so the hypothesis is accepted on the .05 level. Another way of carrying out 
procedure 7 involves computing 


aot ME Ly , Where c=1+ rt 


is distributed as x°(4) (see formula (17) in procedure 7 and the discussion in 
Appendix 3). Since r = 10, k = 5, we get c = 1 + 1/50 = 51/50 and 


he Ly (2.30) = 2.26. 


and where aati tneh, 


If we wish to have the Type I error = .05, the hypothesis of homogeneity 
will be rejected if 


the upper 5% point of the x’(4) distribution. But the observed value of 
—2rk log L, os 
> 2.26. 
Hence the hypothesis is accepted. An alternative test for homogeneity is provided 
by procedure 9. This involves computing 


U = max T ,/min T; e 
In this case U = 12.99/7.25 = 1.79. 
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We enter Hartley’s table 3 (p. 310) for k = 5, vy = 20. The upper 5% point 
for U is 3.49 (i.e., if we reject the hypothesis of homogeneity for U > 3.49 
and accept if U < 3.49, our Type I error = .05). But we observed U = 1.79. 
Hence the hypothesis of homogeneity is accepted. 


(12) Five items are observed to fail on a life test. The total lives between 
successive failures (first total life is measured from ¢ = 0, when the life test 
started) are 627, 898, 830, 50, and 1028. Use procedures 8 or 9 and Appendix 2 
to test whether 0 is constant throughout. Use the 5% level of significance in 
making computations. 

Solution: Using formulae in Appendix 2, one gets 


log L, = 1/5[log 627 + log 898 + log 830 + log 50 + log 1028] — log (#283) 
= 6.1620 — 6.5317 = —.3697. 
- LD, = 691. 


Reference to the L, tables for n = 3 and k = 5 indicates that the hypothesis 
that @ is a constant should be accepted if L, < .328, if we wish to maintain a 
_ Type I error = .05. The value actually obtained is .691 and so the hypothesis 
is accepted at the .05 level of significance. Another way of carrying out procedure 
8 involves computing 


Since k = 5, we get c = 6/5 and 


—# B Ly > (3.697) = 3.081. 


If we wish to have the Type I error = .05, the hypothesis of homogeneity 
will be rejected if 


Eas " Lis 9.49 


(since (—2k log L,)/c is distributed as x*(4) and the upper 5% point of x’(4) 
is 9.49). But the observed value of 


—2 log Ly 2 Li _ 38. 


Hence the hypothesis that @ is constant throughout is accepted. 
Hartley’s test based on the maximum F distribution (procedure 9) gives 
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Since Pr(U > 202) = .05, we are led to accept the hypothesis that @ is constant 
throughout. 


(13) A life test is carried out and one is given the total lives associated with 
the first twenty failure times. The total lives T(r), --- , T'(720) are 


16 77 131 
18 95 144 
23 120 157 
26 121 195 
27 129 198 


Test the hypothesis that @ is constant in (0, 729) against the alternative that 
it varies in (0, 729). Use procedure 7 (Bartlett’s test), by splitting the 20 observa- 
tions into 4 groups of 5 observations each. 

Solution: If the underlying process is Poisson, then 


2ri{ log 2 2) a i ee TF 4 ME Te ve) + > + eT. — ranal} 


k+1 
6rk 


is approximately distributed as xi_, . In our example the number of groups 
is k = 4 and the number of items in a group is r = 5. Therefore, we have 


2(5)4){ log 198 _ + flog 27 + log (74 — 27) 
BB san ae a enim 


4+1 
1+ FO@ 


+ log (129 — 74) + log (198 — 129)} 


4+1 

1 niaeieaig 

+ 6@@ 

From the x’ table we have that P{x3 < 7.82} = .95. Hence we accept the 
hypothesis that @ is constant at the .05 significance level. 


i 


= 2.20 


(14) Five items are observed to fail on a life test. The total lives between 
successive failures (first total life is measured from ¢ = 0, when the life test 
started) are 68, 1516, 49, 22, and 358. Use Fisher’s test (procedure 10) to see 
whether 7(72 — 7,;) = 1516 is abnormally long. Do this at the 5% and 1% 
significance levels. 

Solution: According to (20), we compute 


y = max (Tn), T(r = 1), +++ TGs = 1) 
T(r) 
In this problem 
Z = 1516/2013 = .753. 


According to a table given by Fisher: 
Pr(Z > .684) = .05 when the number of failures is n = 5. 
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But we observed Z = .753. Hence there is evidence at least on the .05 level 
that T(r, — 71) = 1516 is abnormally long. It is interesting to note that 
Pr(Z > .789) = .01, hence we would not be able to assert that T(r. — 1;) is 
abnormally long on the .01 level. The chances are better than one in a hundred 
of observing a 7'(r. — 7,) longer than the one actually obtained. 


(15) One observes a life test (stochastic process) until 7,, , the time when 
the 21st failure (event) occurs. The observed total lives are: 


T(r;) = 18, T(r.) = 23, Tras) = 44, T(r) = 75, T(rs) = 81, 
T(r.) = 83, T(r2) = 122, T(rs) = 164, T(r.) = 166, T(r1.) = 182, 
T(r1:) = 218, T(ri2) = 247, T(r,3) = 289, T(ris) = 334, 

T(t1s) = 361, T(r) = 373, T(r2) = 377, T(t) = 426, 

T(t19.) = 447, T(t) = 451, T(r21) = 526. 


Are these observations consistent with the hypothesis that the underlying process 
is Poisson with a constant rate parameter? Carry out a two-sided Kolmogorov- 
Smirnov test (procedure 11) at the .05 significance level. 

Solution: Graphical Method (see Figure 3): If the underlying process is Poisson 
with constant rate parameter, then the T(7,;), 7 = 1, 2, --- , 20 are uniformly 
distributed over the interval [0, 7'(7.,)], ie., over [0, 526]. The theoretical c.d.f. 
is given by the equation 


F(T) = 0, z$0@ 
T/526, 0<T < 526 
= 1, T > 526. 


The theoretical c.d.f. is represented graphically in Figure 3 by the line L, . 
The sample c.d.f. associated with the data is given by F,(7'), which is defined 
as follows: 


F(T) =0, for T<0 


are > for T(7;-1) < T(r) < T(r;), j=1,2,--- ,2 
=1, for T(r) > T(r.) = 451. 


Put simply, F,(7’) is a step function, with jumps of 1/20 occurring at the total 
lives T'(7;), 7 = 1, 2, --- , 20 associated with the failure times 7; . 

To test the hypothesis that the underlying process is Poisson with constant 
rate parameter, we wish to see whether F,(7’) differs ‘too much” from F(T). 
Using the table of the Kolmogorov statistic for the 95% point when the sample 
size is n = 20 (see [10]), we can state that if the null hypothesis is true, then 
there is probability .95 that the maximum deviation between F,(7’) and F(T) 
is in absolute value <.2939. Graphically this means the following: 

Draw two lines, L, and L; , parallel to L,.such that the vertical distance 
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PROBABILITY 
* . 


TOTAL LIFE 


Fiagure 3—The Kolmogorov-Smirnov test carried out graphically for the data‘iin Example 15. 


between them and L, is +.2939. Rejection of the null hypothesis at the .05 
significance level will occur if the step function F,,(T) crosses either L, or L; . 
If F,(7) does not cross either of the boundary lines, we will accept. In this 
problem F,(7') does not cross either L, or L; and hence we accept the hypothesis 
that the underlying process is Poisson with a constant rate parameter. 


ANALYTICAL METHOD FOR PERFORMING THE KOLMOGOROV-SMIRNOV TEST 


An equivalent analytical procedure, which is sometimes more convenient, 
is to compute 


D, = lw. |F(T) — F,(T)| 


directly. In our particular problem, this becomes 


Do = max | max {| 556 ~ (Lp) Fis? - (Fn") 


It is readily calculated that D.. = .154.* But from the table of the Kolomogorov 
statistic [10], we have that Pr(D.. < .2939) = .95. Therefore, we accept the 
hypothesis that the underlying process is Poisson with constant parameter. 

Remark: Still another possible solution (see Remark 2 following procedure 8) 
is based on the fact that 


—2 . log Eel = = p> log T(7;) — r log x) | 


* To compute Deo , we readily verify that the maximum deviation is attained at the 10th 
observation. At time T(710) = 182, F,(7) — F(T) = .6 — 182/526 = .5 — .346 = .154. 
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is distributed as x3,_, , when the assumption is made that the underlying process 
is Poisson. In the present problem it can be computed that 


21 
-2| > log T(r;) — 21 log 526 | = 48.18. 
i=1 


From the table of x’ we have that 


Pr(24.4 < xin < 59.3) = .95 


and hence if we use the .05 significance level, we accept the hypothesis that 
the rate parameter is constant, when it is tested against the two-sided alter- 
native that it is not constant. If a one-sided test had been run, we would still 
have accepted the hypothesis that the Poisson rate parameter is constant, 
when it is tested against the alternative that the Poisson rate is decreasing 
in time. This is so because the upper 5% point of xi, is 55.76 and we observed 
48.18. 
Although the null hypothesis is accepted by both tests at the 5% significance 

level, the procedure based on 

> ed | 

X08 | (7, 
is more sensitive. That this should be the case is made clear from considerations 
in [4] and [5]. 


(16) A non-replacement life test is terminated at a preassigned total life of 
4100 hours. The total lives T(7,), --- , T’\(rs0) are: 


151 423 1602 2611 3581 3883 4001 
169 587 1674 2701 3733 3884 4024 
194 766 1722 3084 3738 3907 4040 
195 790 1822 3089 3740 «3914 4043 
203 1131 1989 3192 3776 3929 4069 
205 1138 2009 3326 3807 3933 4070 
218 1162 2234 3351 3812 3946 4070 
231 1306 2438 3428 3814 3957 4071 
269 1493 2458 3496 3829 3977 4088 
148 272 1546 2548 3574 3848 3981 4100 


Test the hypothesis that the underlying process is Poisson with parameter @ 
against the alternative that the parameter 6 has changed. Use procedures 3, 
7, or 9, all at the .05 significance level. 

Solution: In procedure 3 we use the fact that if the Poisson assumption holds, 
then the total lives {7'(7,), --- , T'(rs0)} are uniformly distributed in [0,4100]. 
Because the sample is fairly large, it is convenient to apply the x” goodness 
of fit test, where we are testing for deviations from the null hypothesis that obser- 
vations are uniformly distributed over the interval [0,4100]. Dividing the interval 
[0,4100] into 10 equal intervals we have 
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O- | 410— | 820— | 1230— | 1640— | 2050- | 2460- | 2870- | 3280- | 3690- 


number (e;) 


Observed 
number (0;) 


The calculated value of x’ is 


10 2 
ee 2 (0; = é) _ 

From the x’ table we have that Pr[xf > 16.9] = .05. Hence we reject the 
null hypothesis at the .05 significance level. In procedure 7, we have a choice 
as to how we should subdivide our data into groups. For convenience we choose 
k = 4 groups containing r = 20 observations each. That is, we are interested 
in the total lives associated with the time intervals [0, r20], [720 , Tao], [T40 5 Tool, 
[760 , Tao]. The associated total lives are: T'(r2.) = 272, T(t40 — 720) = 2548 — 
272, = 2276, T (re. — T1) = 3848 — 2548 = 1300, and T (reo — rea) = 4100 — 
3848 = 252. Substituting into (16), using k = 4, r = 20, we get 


2(20) (){og( #100) 1 [log 272-+log 2276+log 1300-+log 252}} 
BB i 


4+1 
1 + a0) 


= 67.07. 


From the x’ table we have that Pr[x} > 7.81] = .05. Hence we reject the 
hypothesis that @ is constant. In procedure 9, we apply Hartley’s U test, i.e., 
we compute 


a max {T (Too), T (To _ T20) , T (Teo — Ta0) 5 T(ts0 —_ Ts0) } 


min {T(720), T(t40 — T20), T'(t60 — Tao), T'(780 —, Too) } 
_ 2276 


= 359 = 9.03. 


Entering Hartley’s table for k = 4, v = 40, we see that the upper 5% point 
of U is about 2.35. (The table does not give the upper 5% point for k = 4, 
v = 40. It does give 2.59 as the upper 5% point for k = 4, vy = 60.) Hence the 
hypothesis that 6 is constant is rejected. 

Remark: The deviations from the null hypothesis are so extreme in this problem 
that one scarcely needs a statistical test. 


(17) A life test is terminated at a preassigned total life of 1000 hours. The 
total lives T'(7,), «++ , (ra) are: 


68 151 205 315 412 625 832 
109 169 218 326 455 700 891 
114 194 231 327 480 704 948 
133 195 269 327 555 757 957 
148 203 272 399 574 = 810 985 
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Test the hypothesis that the process is Poisson with constant rate parameter 
against the alternative that the value of @ increases during the test. Try pro- 
cedures 3, 7, and 11. Test at the .05 significance level. 

Solution: Procedure 3 is based on the result that if the Poisson assumption 
holds, then the total life is uniformly distributed in (0,1000). Since the sample 
is fairly large, the x” goodness of fit criterion may be used to test uniformity. 
Dividing the interval (0,1000) into 5 equal intervals we have: 





| 
| 0-2 200-400 400-600 600-800 00-10 


Expected number (e; ) ae 8 


Observed number (0; ) on oe 


pon 


i=1 


is distributed as x{ . In this problem 
- (0; ie e; P 74 
>! ) 


i=] eC; — 8 
Since Pr{xi > 9.49} = .05, we just barely accept the hypothesis that @ is 
constant over the total life interval [0,1000]. 

Remark: When we applied this test, we did not use the fact that, for the alter- 
native stated, we were testing against an increase in the parameter. This sug- 
gests that an appropriate test might be to compute a rank correlation coefficient 
for the T(r; — r;-:),% = 1,2, --- , 40. 

A test which is more sensitive in this situation is Bartlett’s test (procedure 7). 


From previous considerations we have that if the underlying process is Poisson, 
then 


= 9.25. 


2rk {log Ted) _ } “a {log T(7,,) + + log T'(7;, — aia 


j=2 


k+1 
6rk 


B= 
1+ 
is distributed as x;-_, . 
In our example if we let the number of items in a group be r = 8, and the 
number of groups be k = 5, we have that 


2(8)6){10 og (288) _ 1 ~ | fog 114 + log (205 — 114) 
I sth nba ace alias Michie iarliinligclhndiiiins 


iT @ TE 


+ log (327 — 205) + log (700 — 327) + log (985 — 700) 


1 + aera) 
= 12.346. 


THB 
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Since we are testing against the alternative that @ increases, we have a one- 
sided test. We reject the hypothesis that 6 is constant because from the x’ table 
we have that 


P{xi > 9.49} = .05. 


Finally, we apply the Kolmogorov-Smirnov test (procedure 11). Since the 
life test is terminated at a preassigned total life of 1000 hours, the 40 total lives 
{T(r,)}, 7 = 1, 2, --- , 40 which occurred prior to total life = 1000 hours are 
uniformly distributed over the interval [0,1000], if the underlying process is 
Poisson with constant mean life 6. The theoretical c.d.f. is defined as: 


F(T) = 0, T <0 
= 7/100, O<T < 1000 
= 1, T > 1000. 


The theoretical c.d.f. is represented graphically in Figure 4 by the line LZ, . 
The sample c.d.f. associated with the data is given by F,(7), which is defined 
as follows: 


F(T) = 0, for T <0 


et T(r:-1) < T(r) < T(r), j=1,2,--- ,40 


T(r) > T(t) = 985. 


E 
a 
Z 
gz. 


TOTAL LIFE 


Figure 4—The one-sided Kolmogorov-Smirnov test carried out graphically for the data 
in Example 17. 
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Put simply, F,(T) is a step function, with jumps of 1/40 occurring at the total 
lives & .7;),j7 = 1,2, --+ , 40. 


Since we are testing the null hypothesis that @ is constant against the one- 
sided alternative that @ increases in time, we will use a one-sided Kolmogorov- 
Smirnov type test. From table 1 in [9] we know that 


Pr{ max (F,(t) — F(é)) > .189} = .05, 
0<7T<1000 


when the sample size n = 40. Graphically this means the following: draw a 
line LZ, , which is parallel to Z, and is a vertical distance .189 above L, . Re- 
jection of the null hypothesis at the .05 level will occur if the step function 
F(T) crosses Lz (for, if this happens, maxo<r<iooo (F,(7) — F(T')) does exceed 
.189). In this example F,,(7) crosses Lz and hence we reject the hypothesis 
that the mean life @ is constant throughout the interval [0,1000]. Actually, 


max, (F(Z) — F(T)} = max ib * Tiegh 


1si<40 


24 327 
= 1 _- a = .600 — 327 = .273 
which exceeds .189. As a matter of fact, the 1% point of maxo<r<io00 (Fao(T) — 
F(T)) is .285 and so the null hypothesis would be rejected even at the .01 level. 
For references and descriptions of the theory underlying the above examples, 
the reader is referred to Part I of this paper which appears in the previous issue 
of this journal: Technometrics, Vol. 2, No. 1, February 1960, p. 83-101. 
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Partial Duplication of Response Surface Designs 


O. Dyxstra, Jr. 


General Foods Research Center 


Continuing the author’s earlier work [6] a method is described which requires that 
certain experimental runs of a central composite, second-order, response surface design 
be repeated, thereby providing a more general estimate of the experimental error, at 
the same time providing more reliable estimates of the effects. 

The partial duplication of the factorial portion as well as the partial duplication of 
the star portion has been considered and described. The response surface designs with 
the star portion duplicated seem to have more potential than the designs with their 
factorial portions duplicated or partially duplicated. 


i. Toe ARGUMENT FOR PARTIAL DUPLICATION IN 
RESPONSE SURFACE DESIGNS 


Many, perhaps most, experiments are based on prior knowledge, and what 
is desired is knowledge about the changes in response which occur when the 
present process is changed from some “normal” set of levels. In these cases 
the experimental design is centered according to the experimenter’s best guess. 
Using the presently available second order designs, this “center position” is 
the condition at which the n, replicates are usually obtained. It is frequently 
feared that the variability may increase away from the center point, the resulting 
estimate of error may be too small for a correct assessment of the coefficients 
of the second order polynomial. Thus, since we often have no assurance that 
the experimental error is constant over the entire experimental area, it would 
be sound experimental strategy to obtain duplicates over the experimental area. 
Furthermore, such a study might reveal the existence of heterogeneous error. 

Partial duplication requires additional runs, so that we also obtain more 
precision in the estimates of the coefficients, more degrees of freedom for an 
estimate of the experimental error, and a more powerful test of the adequacy 
of the second order model. 


2. RESPONSE SurRFACE Desians (SECOND-ORDER) 


Before entering into a discussion of partial duplication of second-order re- 
sponse surface designs, it is necessary to review certain aspects of the work 
of Box, Hunter, et al [1, 2, 3, 4, and 5]. 

The number. of variables is denoted by k. For the case of central composite 
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designs (to which attention in this paper is restricted), the total number of 
runs is N = n, + n, + M%o , where 


= 2 is the number of points in a square, cube, hypercube, or fraction 

of the hypercube. These points have coodinates (+1, +1, --- , +1). 

For k > 5 we may take n, = 2°’, and for k > 8 we may take n, = 2°”. 

= 2k is the number of points in a “star” design. These points have co- 

ordinates (ta, 0, --- , 0), (0, ta, 0, --- , 0), --- , , --- , 0, +a). 

m, is the number of points at the center of the design, i.e., with coordinates 
(0, --- , 0). 


With a chosen equal to n? the design becomes “rotatable”, that is, the 
variance of a prediction, V(#), is constant for all points the same distance, p, 
from the center of the design. 

For the central composite designs each of the k variables are at five levels: 
0, +1, and +a. These settings are, of course, a convenient coding of the original 
variables. The np, points at the center of the design provide an estimate of the 
experimental error s*° with n) — 1 degrees of freedom, which can be used for 
statistical tests of significance. These center points serve another important 
function, since they also greatly affect the value of V(g) when p isnot too large. 
Box and Hunter [2] give values of mn) which force V(g) for p = 0 to equal V(#) 
when p = 1. Such designs are called “uniform information’ designs. 

DeBaun [5] and Box and Hunter [2] give methods for blocking second order 
designs, if blocking is necessary. In order to obtain orthogonal blocking, i.e., 
designs for which the block effects are independent of the estimates of the 
coefficients, it is necessary to allot the n, center points to the blocks so that 
the total number of points in the blocks is proportional to the respective values 
of >> a” for each block. 

The true response 7 is some function of the k variables x, , «+ , 2 ; in general 

= 6 (%, , --- , 2%). It is assumed here that this function can be adequately 


ianodaaned by a polynomial function which is at most quadratic in the area 
of interest, i.e., 


(1) n = Boto + Bit: + +++ + Bite + Buti + +++ + Bust 
+ Br2t%%_ + +++ + Be-1,2%n-1% « 


The value of 2, is identically equal to 1. For a given set of conditions the response, 
y. , has mean y, . In general y = » + e, where e are NID(0, o’). 

For the designs under consideration in this paper the matrix of sums of 
products and cross-products, X’X may be written 
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A @ @ --- © BE. FB +s B OL @. «es (xo) 
oe ate Se eee oe fae (x3) 


B+ 00.0 ++ 000 -- (22) 


BOO =: 0 0 --- 0 (2x,) 
CO BR ws. 0 @ -s §g (23) 
et” ves G6 @ (23) 


(xi) 
(symmetric) tien 


(a3) 


(Ta-17) 
Assuming that the variables take on values 0, +1, and -ra, then 
A= Dxa=N=n+nu+m = D+2k+m 
B= Dizi. = D+ 20” forall i 


C= Dict, = D+ 2a forall i 
D= Dawa. =n. forall i+ j 


1,j=1,--- ,k; u=1,---,N 


Although this notation is not used in the previous work, we have found it 
convenient. All cross-product terms involving either x; or 2,2; are zero. 


Box and Hunter [2] have shown that in order to have rotatability it is necessary 
that 


C =3D = D+ 20 


so that a = n? . For a rotatable central composite design the variances and 
covariances of the coefficients are 


V(bo)/o* = (k + 2)D/[(k + 2)AD — kB*] 

V(b,)/o? = 1/B 

V(b;,)/o* = [(& + 1) AD — (& — 1)B’)/2D[(k + 2)AD — kB") 
V(b.i)/o? = 1/D 

Cov (bo , b;:)/o” = —B/[(k + 2)AD — kB’) 

Cov (b;; , b;;)/0* = (B? — AD)/2D{(k + 2)AD — kB’ 
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The variance of a prediction 9 is given by 
Vp) + 2)AD — kB*)/o° 
= (k + 2)D + p(k + 2)(AD — B’)/B + p'[(k + 1)AD — (k — 1)B*]/2D 
where p? = )-*_, 2%. 
The quadratic coefficients will be correlated unless Cov (0;; , 0;;) = 0. In 


order to have an orthogonal design, it therefore is necessary that AD = B’, 
or that 


(7) no = B/D —n, — 1, 


In determining the value of m) which would yield uniform information, Box 
and Hunter coded the variables so that es zi, = N;i = 1, --- , k. In these 
coded units the value of N which equates V(g) at p = 1 to V(g) at p = Ois 
determined. To satisfy this requirement we multiply each z;, by (N/B)* = c. 
In these coded units we have >.™_, 27,23, = AN, i < j, so that c‘D = XN and 
N = B?)/D. Substituting these values into (6) we find that uniform information 
is obtained when 


(8) X= [kK +3 + V9k? + 14k — 7]/4(k + 2) 


Values of \ are given by Box and Hunter [2]. The number of center points needed 
to provide uniform information is 


(9) nm = B/D —n. —% 


In terms of the original coding (0, +1, and +a), V(g) at p = 0 equals V(¥#) 
atp =. 


(6) 


3. PARTIALLY DUPLICATED SECOND ORDER RESPONSE SURFACE DESIGNS 


In this paper eight types of partially duplicated second order response surface 


TABLE 1. Description of Designs 


II 
III 
IV 
V 
VI 
VII 
VIII 


Ne 


Qk 

Qk-1 

DQk-2 

Qk a Dk 

Qk 4+ Qk 
Qk-1 4. Qk-1 
Qk-1 4. Qk-2 
Qk-2 4. Qk-2 
Qk 

Qk-1 

Qr-2 


Descriptions 


Full replication 

Half replication 

Quarter replication 

Factorial duplicated 

Full replication plus half rep. 
Half replication, duplicated 
Half plus quarter replication 
Quarter replication, duplicated 
Full rep., star duplicated 
Half rep., star duplicated 
Quarter rep., star duplicated 


For all three classes of designs D = n,. For classes 0 and 1 we have A = D + 2k +n, 
B = D + 20%, andat = D. For class 2 we have A = D + 4k + m, B = D + 402, anda‘ = 
D/2. Tables 2, 3, and 4 summarize some details of these designs. 
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designs are presented. These, in turn, may be broken into two classes of designs 
which, with the usual Box-Hunter designs indicated as class 0, are described 
in Table 1. 


TABLE 2. 


Values of ne, Na, and no for Rotatable Designs With Full Factorials, k = 2, ... ,8 
(O = orthogonal, U = uniform information) 


TABLE 3. 


Values of ne, Na, and no for Rotatable Designs Where Half Replicates Are Possible 
k = 5, ... ,8(O = orthogonal, U = uniform information) 


Ne 


Ne 
nO) 
ng U) 


Ne 


Na 
nO) 
no U) 


Ne 


Na 
nO) 
no U) 


Ne 
Na 
nO) 
no U) 


a 
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TABLE 4, 


Values of ne, Na, and no for Rotatable Designs Where Quarter Replicates are Possible 
k = 8(O = orthogonal, U = uniform information) 


Type k=8 


IV Ne 192 
Na 16 
nO) 43 
no U) 25 


Ne 128 
Ne 16 
nO) 33 
no U) 20 


Ne 64 
Ne 32 
n&O) 21 
no U) 13 


Ne 64 
Na 16 
nO) 20 
no U) 13 


4. BLocKING 


As in the designs given by Box and Hunter, not all of the new designs can 
be orthogonally blocked and still retain rotatability. If it is desired to have 
orthogonal blocking, then a value of a can be found with a suitable choice 
of m» in each block such that the designs are nearly rotatable. Descriptions 
of some blocking arrangements are given for all three classes of designs in 
Table 5. 

To achieve orthogonal blocking for classes 0 and 1, it is necessary to distri- 
bute the , center points to the factorial and to the “‘star’” designs so that 


(10) (D + .m20)/(2k + eo) = D/2a* 


where 7.9 + Mao = MN. Since N49 = N — D — 2k — n. , we can readily deter- 
mine No by 


(11) no = N[D/(D +2VD)] — D 


If, with a reasonable value of N, it is not possible to obtain a whole number 
for n.o , then having chosen values for N and 7,0 , the value for a may be deter- 
mined by 


(12) a = (D/2)[((N — D — no)/(D + n0)] 
For class 2 the center points are to be distributed so that 
(13) (D + N.0)/(4k + Nao) 7 D/40" 
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TABLE 5. 


Orthogonal Blocking Arrangements for Some Central Composite Designs 

with Approximately Uniform Information 

k = number of controlled variables 

me = D = number of points in cube (factorial portion) 

neo = number of center points for blocks of cube 

Nq = number of points in star (2k for classes 0 and 1, 
4k for class 2) 

Nao = number of center points for blocks of star 

N =m. + neo + Ma + Nao = total number of points 

a = radius arm of star design 


N Blocks** 


> 
3 
3 
= 
oO 


3. 1. 4+2;44+2 
au “¥, : 2(4 + 2);6 +2 
30 2. 28 + 2);8+2 
54 2. : 4(8 + 2);10+4 
33 2. 16 + 6;10+1 
54 2. : 2(16 + 4);12 +2 
90 2. 88 + 1);14+4 
95 2. 4(16 + 3);16 +3 
90 2. 4(16 + 2);16 +2 
24+1);4+3 
4(4 + 1) or 2(8 + 2);6 + 4 
48 + 2);8+6 
&(8 + 1) or 4(16 + 2);10 +8 
2.667(2.632) 3(16 + 2);10 +6 
2.366(2.378) 2(16 + 4);10 +4 
2.404(2.378) 2(16 + 2);10 +3 
2.828 4(16 + 2);12 + 6 
3.361(3.364) 8(16 + 1);14+ 10 
3.694(3.722) 6(32 + 2); 16 + 13 
3.757(3.722) 6(32 + 2); 16 + 14 
3.361(3.364) 8(16+1);16+8 
1.195(1.189) 4+ 3; 2(4 + 1) 
1.414 2(4 + 3); 2(6 + 1) 
1.683(1.682) 2(8 + 4);16+1 
1 664(1.682) 2(8 + 5); 2(8 + 1) 
000 2(16 + 5); 20 +1 
000 48 + 3); 2(10 + 1) 
690(1. 16 + 12; 2(10 + 0) 
000 2(16 + 8); 2(12 + 0) 
366(2.378) 8(8 + 2); 2(14 + 0) 
412(2.378) 4(16 + 6); 2(16 + 0) 


Ono & YY 


_ 


1-II 
1-III 
1-III 
1-III 
1-III 
1-IV 
1-IV 
.1-V 
2-VI 
2-VI 
2-VI 
2-VI 
2-VI 
2-VI 
2-VII 
2-VII 
2-VII 
2-VIII 


— 
WDWAOWWONNAMKFRAOOWA RN WNW OC 


= 
SCOOCON KE NYKNYNYWRWOAWKRAWAKRWNWH NK PND WD 


28 
32 


ONMHAIKIKTEP RP WNHWWDWDNOUAATAAhWNHWWONOAT PP WH 


2. 
2. 
zs 
2. 
2. 
2. 


* When the value of a for rotatability differs from the value of a for orthogonal blocking, 
the value for rotatability is given in parentheses. 

** Given are the number of cube points in a block, the number of center points for the 
block containing the cube points; the number of points in the star, the number of center points 
for the block with the star points. Thus 2(4 + 2); 6 + 2 indicates three blocks, two of which 
contain 4 cube points and 2 center points and the third of which contains 6 star points and 
2 center points. 
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Since n.o = N — D — 4k — n.o we have 


(14) no = N[(D/(D + 4V D/2)] — D 


and 
(15) a = (D/4)[(N — D — no)/(D + n0)] 


Further blocking is possible for many of the partially duplicated designs. 
For class 2, type VI, k = 4, for example, we haven, = 2* = 16, n, = 2(2k) = 16, 
and n, = 9 for uniform information. The design is virtually rotatable if n.. = 8 
and 4, = 1. Three blocks are possible: 


1. A one-half replicate of the 2* plus 4 center points = 12, 
2. The other half of the 2* plus 4 points = 12, and 
3. The 8 star points, replicated, plus 1 center point = 17. 


Another almost rotatable design is obtained with 12 center points by setting 
Neo = 10 and no = 2. Four blocks are possible here: 


1. A one-half replicate of the 2* plus 5 center points = 13, 
2. The other half of the 2* plus 5 center points = 13, 
3. and 4. The star design plus 1 center point = 9 each. 


5. Discussion 


In this section the designs will be compared in certain important respects, 
namely, (1) biasses in the coefficients due to the presence of third order terms, 
(2) power of the test of the adequacy of the second order model (Lack of Fit 
Mean Square/Error Mean Square), and (8) information curves. In addition, 
the problem of testing for heterogeneous error is discussed. The following section 
will discuss some experimental strategies. 

In order to compare the various designs, the designs must first be made 
equal in “size”. Two designs will be considered equal in size if >>, x3, = N 
for each design. This amounts to multiplying each x by ¢c = (N/B)', if the 
values of x are 0, +1, and +a. In coded units the mean squared distance of the 
experimental points from the center of the design is k. The moments as originally 
defined, D = n, = [ttjj] for all i ¥ j and B = [i] = D + 2a’ or D + 4a’ for 
all 7, will be used in the following sub-sections. 


a. Bias Due to Third Order Terms 


The method described by Box and Hunter [2, page 198] is used. For all central 
composite designs covered in this paper it is found that the bias in the linear 
coefficients is given by 


(16) &(b,) = 6B; + (D/BY( Buu + d Bui). 


Equation (16) shows that, for any n, and n, , the bias in the linear coefficients 
increases as 7 increases. If the term ND/B’ (= X for uniform information 
designs) is constant for the designs possible for a given k, then the biasses 
in the linear coefficients are equal for these designs. 
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For design types A, I, and VI none of the coefficients are biassed by the 
presence of @,;; terms, where h # 7 ¥ j;h, 2, j, = 1, --- , k. For k > 6 this 
also is true for design types B, II, III, and VII. For k = 5, however, design 
types B, II, III, and VII include a fractional factorial which has each of the two- 
factor interactions (8;; terms) confounded with a three-factor interaction 
(Bxtm terms). Thus for k = 5, 


(17) &(b:;) = Bi; + Brim ; 
where z = (16/D) (N/B)' for design types B, II, and VII and z = (32/D) (N/B)? 


for design type III. The bias is considerably less for design type II than for 
the other three design types (2 = 0.36 vs. slightly larger than 1.). 


b. Testing for Lack of Fit 


The adequacy of the second order model is tested by comparing the Lack 
of Fit Mean Square with the Error Mean Square. For the types of central 
composite designs given in the paper the degrees of freedom for these mean 
squares are given in Table 6. With ND/B’ constant as in the case of uniform 


TABLE 6. Degrees of Freedom for Lack of Fit and for Error 


Type of Degrees of Freedom Degrees of Freedom 
Design for Lack of Fit for Error 


I 2-+n—1 
II 2* — k(k — 1)/2 214) — 1 
VI 2k +n — 1 
A : No — 1 


III 2-1 + no — 1 

VII 2-1 — k(k — 1)/2 2k + no — 1 

B no 1 

V 2*2 — k(k — 1)/2 2*27 + —1 
2k +n — 1 

Cc no — 1 


IV [depends on the design] 


information designs, for a given k the new designs, compared with their Box 
and Hunter counterparts, have more degrees of freedom for error, so that the 
test for the adequacy of the second order model will have more power. 


ce. Information Curves 
The information, (9), is defined as 
(18) I(g) = °° /V@), 
and the information at the center of the design, J,(9), is given by 
(19) (9) = 0° /V(bo) = ro +2. +2. — kB’/(k + 2) D 
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Certain values of I,(g) are given in Tables 7a and 7b. It is seen in these tables 
that only in a few designs the external points contribute substantially to the 
information at the center of the design. 


TABLE 7a. Io (9) for Design Types A, I, and VI for k = 2, 3, and 4 


Design k 
Types 3 


No No ate 0.01 no 
no + 0.34 no + 0.40 no + 0.92 
no + 0.34 no + 0.80 No + 0.92 


TABLE 7b. Io (9) for Design Types I, II, III, A, B, VI, and VII for k = 5,6, and 7 


Design 
5 


No + 2.57 
no aod 1.06 
no ae 0.12 
no + 0.29 
No + 0.57 
VII no + 2.70 


Information curves are given by 


i 2n[(k + 2)d — KIN 
20) IO = FEF ON F BEF DA — De + (KE FDA—-E&—DIo 


where \ = ND/B’ and there are N observations. Since the value of \ depends 
on k, I(g) is a function of k, N, and p. Comparing two uniform information 
designs with N, and N, runs, the information curves for every p will be in 
proportion N,/N,.. Thus it is seen that any of the new designs, compared 
with the Box-Hunter designs, will give increased precision in direct proportion 
to the number of experimental runs. 


d. Testing for Heterogeneous Error 


In the new designs duplicates are obtained at the center of the design and 
either in the hypercube or the star design. Thus two estimates of error may 
be obtained, which we hypothesize are equal. The agreement between the 
two estimates may be tested by a two-tail F test. 

It would seem that another method may be used if the duplicates are obtained 
in the hypercube. A factorial analysis of the absolute values of the differences 
between the pairs of duplicates should reveal the presence of heterogeneous 
error, if such exists. 
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6. EXPERIMENTAL STRATEGIES 


In experimental situations where the experimental error is fairly small, the 
designs of Class 2 will be found useful, since they permit partial duplication 
at minimum cost. 

When little or nothing is known about the magnitudes of the effects relative 
to the experimental error, the designs of Class 1 might be used. A logical pro- 
cedure would be to follow a sequential method of experimentation as suggested 
by C. Daniel [7], starting with a relatively small fractional factorial, except 
that center points are included in each block of runs. If no effects are found 
the experiment may be abandonded for a new experiment with the levels spread 
further, or additional blocks are run. Running additional blocks is continued 
at least until the factorial has two-factor interactions confounded with higher 
order interactions. 

For k = 5, for example, we could begin with a one-half replication of the 
factorial, including some center points. Depending on the results (1) the star 
design and some center points can be run, to complete a type B design, (2) 
the block can be repeated and the star design and center points can be run, 
to complete a type III design, (3) the other half of the factorial can be run, 
again with some center points. If, after completing the factorial, still more 
information is needed, all or half of the factorial may be duplicated, again 
with center points. Eventually the star design will be run to complete the 
central composite design. The result will be a type I, type II, or type A design. 
In general we would prefer to use the Class 1 designs for this purpose, because 
the Class 2 designs require too many center points for orthogonal blocking, 
and also because, having started with the factorial, the radius arm of the star 
design for the Class 1 designs will have to be larger in order to achieve a rotatable 
composite design. 
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A Rank Sum Test for Comparing All Pairs 


of Treatments* 


Rosert G. D. STEEL 
Mathematics Research Center, United States Army, 


Madison, Wisconsin 


A multiple comparison rank sum test, for the simultaneous comparison of all pairs 
of treatments in a one-way classification with equal numbers of observations, is pre- 
sented. An example is worked and tables of critical values are given. Computation of 
probabilities for the general case of unequal numbers of observations is considered and 
means, variances, and covariances are given for this case. 


1. INTRODUCTION 


Recent research in multiple comparison procedures has been concerned with 
non-orthogonal sets of comparisons. Certainly some of the enthusiasm with 
which such problems have been approached is due to the natural desire to 
detect real differences among treatments even when this involves testing effects 
suggested by the data, for example, testing the difference between the greatest 
and least average treatment responses. Among such multiple comparison pro- 
cedures, that which, in the framework of the analysis of variance, tests all 
possible pairs of treatment means, has received most attention. 

In connection with multiple comparison procedures, several error rates have 
been defined. Two of these are the per comparison error rate and the experi- 
mentwise error rate. The former is defined as the ratio of the number of false 
conclusions drawn when the null hypothesis is true to the total number of 
conclusions drawn. This is the error rate customarily used when comparisons 
are independent, for example, in factorial experiments. An experimentwise error 
rate is defined as the proportion of experiments in which at least one false 
rejection of the null hypothesis is made. 

Criticism has been levelled at both error rates. For a per comparison error 
rate, the probability of at least one false rejection of the null hypothesis in- 
creases with the number of treatments in the experiment. For an experiment- 
wise error rate, the probability of detecting a real difference decreases as the 
number of treatments in the experiment increases. How applicable these facts 
are as criticisms depends upon such considerations as the size of meaningful 
differences, the importance of detecting them or failing to detect them or falsely 


detecting them, the number of treatments, and the choice of significance level 
for the chosen error rate. 


*Sponsored by the United States Army under Contract No. DA-11-022-ORD-2059. 
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This paper presents a fixed rank sum test for comparing all pairs of treatments 
in a completely random design. Error rate is experimentwise so that the test 
is an analogue of Tukey’s (3) fixed range test; this test is described by Federer 
(1) and ‘others. For a fixed rank sum test, it is necessary to compute the prob- 
ability distribution of a certain minimum rank sum. One can then test all pairs 
of treatments as in Section 3. The computation of probabilities is discussed 
in Section 4 and several distributions are presented in Tables 1 and 2. Means, 
variances and covariances are given for the general case. 

Kruskal and Wallis (2) have given a rank test for one-way classifications 
which is an analogue of Snedecor’s F-tests. Thus if it is concluded that the 
test criterion cannot reasonably be attributed to chance and the null hypothesis 
of no differences, one would probably conclude that at least the most extreme 
difference was evidence of a real difference between the corresponding popula- 
tion means. One could hardly conclude more without resorting to further testing, 
thus using a result-guided test procedure. The rank sum test presented here 
permits comparison of all treatment pairs without recourse to a result-guided 
procedure. 


2. PROCEDURE 


LetiX, ,.--: , X, be continuous random variables measuring a single charac- 
teristic on each of k treatments. Represent their cumulative distribution func- 
tions by F, , --- , F, . We wish to test the null hypothesis 


H,:F, = eee = F, 


using n observations from each population, against the alterrative 


H,:F; A F;, at least one?, j pair. 


(More generally, one might take n, , --- , nm observations. However no tables 
have, as yet, been prepared for such a test.) 
To test Hy , 


1. Jointly rank the X,’s and X,’s, X,’s and X;,’s, --- , X,’s and X,’s, X,’s 
and X,’s, --- , X,’s and X;,’s, --- , and finally the X,_,’s and X;,’s, assigning 
rank 1 to the least observation and 2n to the greatest in each ranking. 

2. Add the ranks of the X,’s when ranked with succeeding X,’s to give rank 
sums 7}. , --* , 7, ; add the ranks of the X,’s when ranked with succeeding 
X,’s to give T.,, --: , Tx 5 --* ; finally T,-_, , is obtained. Compute conjugates 
T!; = (2n + 1)n — T,; . The smaller of 7;; and 74; is to be used. 

3. From Table 2 and using the predetermined significance level for k and n, 
extract the value of the significant rank sum. (For a one-tailed test, Table 1 
gives significant rank sums. Such a test would be applied only if one were pre- 
pared to rank all treatments in stating H, ; if the treatments are then numbered 
according to their rank under H, , it is only necessary to compute the 7',,’s.) 

4. Compare the tabulated rank sum with each of the computed rank sums, 
declaring any computed sum that is as small as or smaller than the tabulated 
sum to be significant. Error rate is experimentwise. 
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It may be noted that the procedure can be carried out by ordering the com- 
plete set of observations from smallest to largest and then ranking as described 
in step 1. 

Let us investigate the effect of heterogeneity of variance in the parent popula- 
tion on the sensitivity of the proposed rank sum test as a test of location. Consider 
the extreme situation of two populations so different in variability that observa- 
tions on one of them never lie within the range of observations on the other. In 
other words, assume one population to have variance which is virtually zero and 
the other to have a sizable variance. Then for n observations in each sample, the 
only possible values of the rank sum criterion are n(n + 1)/2 by steps of n 
to n(n + 1)/2; the probabilities are given by the binomial distribution with 
p = } and exponent n. 

The most extreme sequence in this extreme situation occurs when all observa- 
tions from the more variable population fall on the left (or right) of the observa- 
tions from the less variable population. Let us compare the probability of this 
event with that for the same event and identical populations. The ratio is 


1 (Qn)! 1 (Qn + eee V 2 
2" (nly? 2 (n+ 1 eae 

_ 1 (2n)"**(1 + 1/2n)"*4e 

ort nr + 1/n)?*) Vr 


using Stirling’s approximation for factorials. Since lim,.. (1 + 1/n)" = e, 
this expression approaches 2”/*/an. Hence the ratio of these probabilities 


almost doubles itself for large n each time n increases by unity. For n = 3, 4, 5, 
exact values of the ratio are 2.5, 4.375 and 7.875 respectively. 

It is unfortunate that the most extreme value of the test criterion is more 
frequent with heterogeneity of variance than with homogeneity. However, we 
have examined a most extreme case of heterogeneity and the most extreme 
term in this situation. In addition, we must bear in mind that means and vari- 
_ ances are positively correlated in many experimental situations and a test 
which, in such a circumstance, detects differences between variances, even though 
they be attributed to differences between means, would be a useful test. 


3. EXAMPLE 


The data used to illustrate the technique are a random sample from a larger 
set. The variable is air permeability of soil in units of square microns and meas- 
ures the rapidity with which air flows through soil; low values are associated 
with compacted soils. All measurements were made on the same day. There 
is some evidence that the log normal distribution may be appropriate. There 
were three soil treatments. 

For these data, the minimum value of the test criterion is 14, the conjugate 
of an observed rank total. From Table 2, we find that values less than or equal 
to 14 or equal to or greater than 22 occur with probability .625. There is no 
evidence of differences among treatments. (This conclusion does not hold for 
the complete set of data, 16 observations on each treatment.) 
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Data obtained through courtesy of Professor Champ B. Tanner, University of Wisconsin. 


For the Kruskal-Wallis (3) procedure, all treatments are compared simul- 
taneously but without benefit of multiple comparisons. Error rate may be 
described as experimentwise. The test criterion has the value H = 1.885, H 
being distributed approximately as x” with k — 1 = 2df when the null hypothesis 
is true. We find .70 > P(H > 1.885) > .60, essentially as found before. An 
important difference in the test criteria is that the former, a multiple comparison 
procedure, will point out which comparisons show significance. 

If the Wilcoxon (4) two-sample test is used, then error rate is per comparison. 
The probability of obtaining a value of the test criterion less than or equal 
to 14, the observed minimum, or equal to or greater than 22, the complement 
of 14, is .342. Comparing this probability with an experimentwise error rate, 
we find that values less than or equal to 12 or equal to or greater than 24 occur 
with probability .260 while values less than or equal to 13 or equal to or greater 
than 23 occur with probability .417. In other words, the per comparison prob- 
ability level of P = .342 for the minimum rank sum of 14, if chosen as an experi- 
mentwise error rate for the proposed multiple comparison test, would call for 
an observed minimum rank sum of 12. However, this value is more severe 
than required, P = .260, while 13 is not sufficiently severe, P = .417. 

Since the only tie occurring in these data is for observations on the same 
treatment, no problem is introduced. When ties occur for observations in dif- 
ferent treatments, it is suggested that they be assigned consecutive ranks in 
such a way as to increase or decrease the value of the test criterion, the decision 
depending upon the seriousness of the consequences of a wrong conclusion. 


4. MATHEMATICAL APPENDIX: DISTRIBUTION OF min T’,; 


Let X, , --: , X, be continuous random variables measuring a single charac- 
teristic in k populations on which are made n, , --- , n observations respectively. 
There are (>, n;)!/(n;!) possible permutations of such a set of observations, 
all permutations being equally likely under the null hypothesis. Hence if all 
possible permutations are considered, if ranks are assigned to the observations 
from the left by pairs of variables, and if the sum of the ranks associated with 
the variable of lower subscript is computed for each pair, then we generate 
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TABLE 1 


(one-tailed—observed ) 
Probability that minimum 77; is less than or equal to tabulated 7;; for k = 3. 


the distribution of (Ti. , «+: , Tix, Tos, -** » Tox, *** » Te-1a)- The first sub- 
script on any 7';; indicates the variable for which ranks are summed. From 
this joint distribution, we can obtain the distribution of min 7’; . 

Let W(Ti2 , «++ , Te-12| m , +++ , M) represent the number of ways of ob- 
taining (T,. , --- , T-1,,.) by summing the ranks of the X,’s, --- , X,-1’s when 
ranked with each other and the X,’s as described in the preceding paragraph. 
If an X, is the last observation in a permutation, then the ranks n, + n,,°-: , 
m + nm, are in the rank sums 7), , --- , 7, respectively and in no others, and 
the remainder of each rank sum, namely 7,; — n, — n; , must be obtained 
by drawing n, — 1 ranks from n, + n; — 1 ranks, j = 2, --- , k. If an X, is 
last in a permutation, then the ranks n, + m3, °:* , M2 + 7 are in the rank 
sums T',, , --- , T'2, respectively and in no others, and the remainder of each 
rank sum, namely 7; — n. — n; , must be obtained by drawing n, — 1 ranks 
from n, + n; — 1 ranks, j = 3, --: , k. Finally, if an X; is last in a permutation, 
no T';; is affected. The number of ways in which (Ti. , «++ , Tx-1.4) may be 
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TaBLy 2 
(two-tailed) 
Probability that either minimum observed 7;; or conjugate of maximum observed 7’,; 
is less than or equal to 7; for k = 3. 


n=l n=2 n=3 


re ee 
1 1.0000 


.0062 
0122 
.0240 
.0413 


obtained is the sum of the numbers of ways in which each of these outcomes 
may occur. Thus we have equation 1. 


WT. , hs » Troi | gas 9) 


= WT — m — mM, °°? > Tu — My — mM, T2334 °° 


ag 4 


+ W(T 1. pee » l'y-2,k eerie — M-1 — % | n, yy eae 1, m) 
+ WT. ,; oh 5 » Tenia | gS see oe 1) 


Equation 1 may be used recursively till any term may be evaluated by in- 
spection. Probabilities are obtained by dividing values given by equation | 
by the number of possible permutations. 
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For k = 4 and n = 1, it is readily demonstrable that equations 2, 3, and 4 
are valid. It is then clear that they hold, where applicable, for all k and n = 1. 


AT ;) = ET.) = n(n; +; + 1)/2 (2) 
o(T;;) = nnn; +n; + 1)/12 (3) 
o(T iT) = mnan;/12 = o(T4;T i) 

o(T 4753) = —nnn,/12 (4) 
o(T 7 ;) = 0 


(Recall that the first of any pair of subscripts on a T' is always the smaller.) 
It is only necessary to show that equations 2, 3, and 4 satisfy the obvious re- 
cursion formulas based on equation 1 in order to prove their general applica- 
bility. For convenience in handling subscripts, equations 2, 3, and 4 will be 
shown to be valid for E(T',2), o°(Ti2),¢ (T2713), (T'12T 23), o(T'137'23) and o(T 27's). 


E(T 2) - zu T12W(T 12 go SPs Tr-1.2 | m en sm)/ De W(T 12 eos 4 Ti-1.% 


[my +++ , ms) 
= {> (Ce —m —m) + (un +2.) WT» —m —m,-+ Te —m — Mm, 
Ta, *** fesslh ~ hey jal 
+ DVTeW Tie, +++ Tu, Tos — te — Ms, *** » Toe — Me — M|, 
Tsu, °°* » Te-r.n | ,M% —1,%,°°* ,M) 
de vee 
+ DET WTrs y +++ 5 Te-a.e y Trai. — Maa — Me | My +++ Ma — 1,m) 


+.Ty.W(Ti2 °° y Ti-ae | m 9 °°* y Men » Ne — ») / and! 
- | + 1m +m) m5 mlm + 1) te 
2 > 1; y > 1; 
M(t + e+ 1) ms Milt + +1) m1 
+ os + > 2 


2 n; 


4. mlm + ms + 1) | (Qin)! _ mim + me + 1) 
2 


2 nN; a(n;!) 


Hence, equation 2 is valid. 

The first term in the above expansion has been considered as two terms 
for purposes of summation. For m, and nz , the set of 7.’s is >.%., 7 by steps 
of one to ae (m + n, — 2) and for n, — 1 and nm, it is Ze zt by steps of 
one to >.%2;' (m + me — i). Hence in summing over the set of T'2’s, which 
is based on n, and nz , the derived set of (T:2 — nm; — 1)’s includes all the 
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values appropriate to the pair (n, — 1, n.), which appears in this first term, 
plus additional values which have zero probability in the distribution arising 
from (n, — 1, 2). 


E(T?.) = >> T2,W(Tia , -** y Te-a.n | y *** yM)/ Dy WT 0? Thea 
| aiees 7) 
= {> [Tis — m — 1)? + 2 + m2)(Tr2 — m1 — M2) 
+ (n, + 2)" JW(T 2 — mm — Me, 0s Tue — Mm — Mm, T35 °°? 5 Tee 
|m, ~ 1,2 ,°** ,m%) 
+ >> T2W(Ti2 ++ Tu, Tos — M2 — M3, °° » Tx — Me —™, 
T sa 5 °°° » Te-a.n | > M2 — 1,73, °° , mm) 
+ >> TWTie » °°» Te-a.e » Te-1,8 — Meni — My | , -°* M1 — 1, m) 


+ x Ti2W(Ti2 poet y Tena | m, 9° °° y Nena» Me — 1} / Rend! 
ae (n, — 1)n.(m. +) _ m es a as 1)*(n, + Ns)” NM 
es 12 = nN; 4 s: nN; 
2 + ne)m (rm — I(r +) , miu + 2)" 
+ Nn; 2 + a nN; 
+ E — 1I)m +n) {. ni(t + tm) (don — 1)In._x(n;}) 
12 a(n; !) (>> n)! 


kK +m +1), nim + + uy") (dons — Ing x(n!) 


12 mn!) (don)! 


Mine(Mi + m2 +1) mer i. ni(n: + m2 + 1)? mer 
12 > 0; ‘te 

4 NyNo(N + 2 + 1) Nt ag ni(m +n, + 1)? nN 
12 ~ ie 4 > 2: 


_ mn(n + m + 1) m ni(n, + m2 + 1)? 
12 4 


Now, 
o°(T12) rr E(T 2) any E°(T'2) 


dn nyn(n, + Ne, + 1) 
12 


and equation 3 is proven. 
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BGT) — Do Pick aW le 3 °° * Tera 
|r, e+ yM)/ >, Wie 5 ** » Tentse |My °° * » Mm) 
= {D0 (Tie — m — m)(Tis — ms — me) + (rs + 2)(Tis — m1 — 1s) 
+ (my + 1s)(Ti2 — m1 — M2) + (rm + 12)(ri + 15)] 
~W(Ti2 — Mm — Me, *** y Te-ae | % — 1,2, *** Mm) 
+ ST uTisW(Ti2 , +++ » Tix, Tos — M2 — M3, *** » Toe — Ma — MH, 
Tse, °°* > Te-a.e | » M2 — 1,3, °°* , Mm) 
+ DTT aW(Ti2 5 ++ Te-o.e y Te-t.2 — Mea — ™ 


My, °°? 5 M1 — 1, m) 


(don)! 
+ 2. T2T'13W(T 12 a Tr-1,% | m4 9 ° °° 9 Mp » hy — »/ at 
_ (i — 1)nnn, 4 (n1 — 1)’n,(n; + ne)(n, + nz) 
“ 12 - nN; + = n; 
4 Qa = V(r + ra)(rr + Ma) 4 (tr — Vals + ma) (nt + 1a) 
2 a nN; 2 z Nn; 
$ ny (Mm + n2)(n1 + Ns) 


nN; 


+ n(n. — 1)nns 4 nyns(m + n)(m + 13 + 1) 
12 z nN; 4 = n; 


+ nyn(nz — 1)ns 4 nins(m, + m2 + 1)(m + ns) 
12 2 Nn; t >» Nn; 


ni(r, +m. + I(m +ns+1) m% 
4 


Nn; 


NyNeNg Ms 


+... 


Mifat _Te ni(m +m + 1)(m +n3+1) m 
13 > nN; 4 a n; 


= Malay ni(n: + m2 + 1)(m, +73 + 1) 
ae 4 


+ 


o(7 273) = E(T 27 13) fo, E(T, 2)E(T;s) 


_ UN2Nz 
12 


which was to be shown. 
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RE. = 5 PTW, -** | Tena 5 8 Sy Wey Tans 
| m yee a) 
{> (Tie — m — M2) + (rr + 2) ]T 23W(Ti2 — ™ — M2, °°? y Trae 
|m, — 1,2, +++ ,m) 
D TrTesW(T 2, ++ » Tix, Tes — le — My °** » To — M2 — M%, 


Ts °°° » Ta-1,0 | »M_ — 1,3 ,°** » M%) 


Do TiTisW(Ti2 y +++ Te-ae |My ++) Mana Me — 1p) / Pana! 


a(n;!) 


— C= Vranats , (rm — Vrinalr, + ra)(M + Ms + D) 
12 ee Nn; 4 2 nN; 
NN (N, + N2)(M2 + 3 + 1) 
2 > '§ nN; 


23 (Ne a 1)nonz 4 (Ne bot 1)n.(n; + Ns)(Ne + Ns) 
12 zs nN; 4 te nN; 


- 


4 NN(N, + No)(N2 + Ns) 
2 ee Nn; 


_ Mitta(tg — Ung , mnsna(rs + m2 + 1)(ma + ns) 


nN; 


_ Musas NMNN(N, + Ne + 1)(ne + 3 + 1) 
12 2 nN; 4 z nN; 


NNNsNe , NNM(N, + Ne + 1)(n2 +73 + 1) 


_ _ MNgNs NN(n, + N + 1)(n2 + 3 + 1) 
satel’ “i 4 


o(T',37'23) = _— 
which was to be shown. 

The two remaining expected values are found similarly and equation 4 is 
seen to be valid. 

Certain symmetry relations hold and are helpful in checking and computing 
tables required in conjunction with the test procedure. Consider a permutation 
of n, X,’s, ++, % X,’s. Ranks are assigned and the test criterion (71. ,--+, Tx-1.s) 
is computed. Corresponding to this permutation is one consisting of the same 





n 4 is 


outing 
tation 
",-1 a) 
. same 
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sequence of X’s in reverse order. When ranks are assigned, the test criterion 
becomes (Tf, , --* , T-1,,) where 77; = n(n; +n; + 1) — T;; and both T;; 
and T{; are the sums of n; ranks, n; < n; . Hence, we have equation 5. 


W(T12 , °° » Ta-1.0 | = >) = W(Ti2 , °° » Ties. | % es » Mx) (5) 


It is clear that for n, = --- = m = n, say, we may interchange subscripts 
and derive additional relations. This process together with application of equa- 
tion 5 lead to equation 6 for k = 3. Clearly, such formulas can be generalized. 


W(T 12 ’ T 13 ’ T 23) = WT? ? Tis ? T%s) 


- W(T 23 ’ Tis ’ Tis) = W(T 23 ’ T 12 ? T 13) 
- W(Tis » Ts ’ T 12) = W(T1s ’ T'23 ’ T{2) 


= W(Ti; ’ T*. ’ Ti.) = W(T23 ’ Tis ’ T 12) 
Bs W(TI2 ’ T 23 ’ T 3) = W(T,2 ’ Ts ? T ts) 


it W(T1s ’ Ti2 , T's) W(Tis ’ Tie ’ T 23) 


In constructing the joint distribution, little use was actually made of equation 1 
except for some checking. Instead, an enumeration procedure and application 
of equations 5 and 6 were found to be more amenable to what was basically 
a hand operation. Tables 1 and 2 were then obtained from the joint distribution 
by simple enumeration. 


(6) 
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The Percentile Points of Distributions 
Having Known Cumulants 


Sm Rona.p A. FisHER 


AND 
E. A. CornisH 


In an earlier study of the uses of moments and cumulants in the specification of sta- 
tistical distributions, the authors developed explicit asymptotic expansions, expressing 
any desired percentile point of such distributions in terms of known cumulants. 

The general formulae are now presented as far as the sixth adjustment, based on the 
eighth cumulant, and also numerical tables showing the coefficients of all terms for ten 
chosen levels of significance over the range 0.5 to 0.0005 (single tail), together with the 
first five Hermite polynomials and tables for the common tests of significance, x?, ¢ and 
z, at the same levels. ' 


I. INTRODUCTORY 


In 1937, in a study of the uses of moments and cumulants in the specification 
of statistical distributions, the authors (1937)* were led to develop explicit 
asymptotic expansions, expressing any desired percentile point of such distribu- 
tions in terms of known cumulants. The general formulae were presented so 
far as the fourth adjustment, based on the sixth cumulant, and also numerical 
tables showing the coefficients of all terms for nine chosen levels of significance 
over the range 0.25 to 0.0005 (single tail), together with the first five Hermite 
polynomials at the same levels. 

As an illustrative example, the cumulants of the z distribution were expressed 
in terms of 1/n, and 1/n, , the reciprocals of the two numbers of degrees of 
freedom, and the rapid convergence at the 5% point exhibited for the case 
n, = 24 and n, = 60. In the intervening period, the formulae have frequently 
been found useful, either for calculations of higher accuracy in the case of 
functions already tabulated, or for values outside their range, or especially 
with tables of multiple entry to supply intermediate values more accurate than 
can be obtained by interpolation (Fisher 1941, Goldberg and Levine 1946), 
or for cases where no tables existed (Johnson and Welch 1939). 


Il. Tue Meruop or Expansion 

The several steps in the method of expansion are set out in Sections 7 and 8 
of our previous paper. Here a brief outline of these steps will be sufficient, 
as the principal purpose of this note is to extend the formulae and tables to 
the sixth corrective term, so widening the range of useful application. 

If the element of frequency in the distribution of a variate & is f(¢)dé, the 

* Attention is drawn to this reference, as it has been repeatedly misquoted. 
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*(-4)}, 


when acting on the frequency function, is to increase the rth cumulant by a,, 
but leaves the distribution otherwise unchanged. This important operationa! 
property of the cumulants is basic to the method, the essential steps of which 
are as follows: 


effect of the operator 


(i) If the cumulants x, , x. , «: + of the distribution of ¢ are expressible in power 
series of the reciprocal of some number n, the frequency element may be re- 
presented as 


Be av”? : q ie ; 3/2 ad gape hy a - = 5/2 £. 


= 700 fae ao es e - genres dé 
V 2rv 
1 


where the coefficients a and c are of order n~”””, b and d of order n~’, e of order 
n°”, f of order n~* --- , and are related, respectively, to x, and xy, Kk, and 
Kk, ks, Ke, °** , and m and» are the mean and variance of a normal distribution 
chosen for convenience. 

(ii) Expanding the operator and integrating, the frequency less than m + év'”” 
may be expressed in terms of the corresponding normal probability integral 
and a series of adjustments of decreasing order of magnitude. 

(iii) If z is the normal deviate at some chosen level of probability, and & 
the corresponding deviate of the distribution under consideration, the difference 
& — x may be found by equating the expression for the probability that the 
variate has a value less than é to 


(1) 


1 : 58/2 1 2/2 2 
vale de — Feetn{e— a) — Le at 
+2¢-2)'@ — 1) -e- ae — 38) + vo} 


(2) 


in which the coefficients are the Hermite polynomials. By considering the terms 
of each order of magnitude in succession, we may develop an expansion for 
& — zx in terms of successive polynomials in é. 

(iv) The expansion for § — x is converted to a much more useful expansion in 
terms of z, the values of which are known in advance, so obtaining the per- 
centile deviate ~ explicitly in terms of the normal deviate. 

The adjustments to the normal deviate having the required probability 
integral are set out in tabular form below. Adjustments V and VI are new with 
this paper. 

Table I gives the numerical values of the first seven Hermite polynomials, 
over the same range of percentiles as given previously, and Table II gives 
the numerical values of the polynomials in the several adjustment terms. 
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III. AtcgEesraic Forms orf THE EXPANSION FOR COMMON 
Trsts or SIGNIFICANCE 


In general applications of the formulae given above, the coefficients a, b, c, - - - 
can be evaluated numerically, and their values substituted, but to facilitate 
the use of the formulae with those statistics in frequent use, the substitution 
has been made algebraically to give the general form of the expansion in terms 
of the parameters of the distributions concerned. 


Adjustment Coefficient Polynomial in x 


1 
2’ —1 


x 
x — 3x 
— (22° — 52) 


—(a2? — 1) 

za — 62° +3 
—(a* — 52° + 2) 
12x* — 532” + 17 


—2£ 

—(2? — 32) 

5(22° — 5x) 

a’ — 102° + 15¢ 

—(2a° — 172° + 212) 

— (32° — 242° + 292) 

142° — 1082° + 1072 
—(252x° — 16882 + 15112) 


2s —1 

—(x* — 62” + 3) 

x‘ — 52° +2 

—(122* — 532” + 17) 

a — 15a* + 452” — 15 

—(x° — 132* + 332? — 9) 

—(a° — 12x* + 292” — 8) 

162° — 1812* + 3932” — 90 

12”° — 1292* + 2712” — 64 

— (802° — 8032* + 1513z” — 304) 
960x° — 8937x* + 15062x” — 2651 





PERCENTILE POINTS OF DISTRIBUTIONS 


Adjust- 
ment Coefficient Polynomial in x 


VI b° x 

bd = B(a* — 32) 

be? §=—s — 35 (22° — 5x) 

bf —(x° — 102° + 152) 

bee = 7(2a” — 172° + 212) 360 

bd = 78a” — 24x° + 292) 768 

b’d —(142° — 1032° + 1072) 64 

bet =—s-:11(2522° — 16882° + 15112) 15552 

h x’ — 2ix° + 1052° — 1052 40320 

cg —(2x" — 372° + 1602° — 1352) 5040 

df —(x" — 172° + 692° — 572) 1152 

é —(2x" — 332° + 1322° — 1082) 3600 
18a” — 2932° + 11002° — 7952 5184 

cde 182” — 2732° + 9742° + 6952 1440 

d 9x” — 131z° + 451° — 3212 3972 

ce —(3962” — 5708a° + 187552* — 118112) 19440 
—(5942" — 81932° + 260062° — 16367z) 13824 

cd = 148” — 670042° + 1952592° — 1095532 62208 

c —(1544402" — 18876842° + 50337142° — 2542637z) 4199040 


(a) x” distribution 
If n is the number of degrees of freedom 
x =n 
+ Vaev2 
+ 2 — 1) 


rd 
i (8 + 142’ — 32 


405 
4 1 (= + 2562° — £334) 
nVn 4860°/2 


as (28 — 2432 — 9232° + 1472) 
n 25515 


n 


pas (27588. + 43532° — 2895172* — aeiii2) 
nvVn 9185400 V2 
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(b) ¢ distribution 
If'n is the number of degrees of freedom 


t= ¢ 


1 (x + *) 
= n ( 4 
en (5 + 162° + az) 
n* 96 
is i, (32 + 192° + 172° — 182) 
n° 384 


eb 1, (a + 7762" + 14822° — 19202° — 2452) 
n* 92160 


ack 1 (2 + 3392° + 9302” — 1782a° — 7652° + 179552), 
n 368640 


The number of terms and the orders of magnitude differ in this expansion 
because it has been derived from the original expansion of Student’s. integral 
in powers of n~’ as given by Fisher (1926). 


(c) 2 distribution 
If n, and n, are the degrees of freedom, and 


ie + 2) 
(ES 
(= — 32x x + us) 
v V5 ne 8e) 4 o =( 72 
ai {3 (« + 92? + 8) if (te! + 72? — is) 
7 120 o 3240 
5 3 
+ BC a) +) a 


+ & (Oe = tae 10180) Lay (2 — 25a — 17727 + a 
o 155520 20160 


4:34 A= + 1012* + 1172? — anh & , (1 + 513a* + 8412? — 2500) 
o 


90720 1632960 











b) 


ion 
ral 


(8c) 










PERCENTILE POINTS OF DISTRIBUTIONS 


i vi{ (z + 72° + 72° + 1058) 
2\° 21504 


+ so'(S0le! + 105112° + 301512° + Sep) 


4838400 
- (sz + 45072° — 829332° — 2643632) 
o 43545600 
_# (ue + 553832° — 3688972° — i21vere)\ 
o 1175731200 ; 


™ As an illustration of the accuracy, we may use the example given previously. 
When n, = 24 and n, = 60, the 5% value of zis 0.26534844, and the asymptotic 
expansion yields the following values: 


Order of ? 

magnitude Successive terms Successive totals Successive errors 
0 0.2809 1224 0.2809 1224 0.0155 6380 
1 — 196 0643 2613 0581 — 40 4263 
2 44 6851 2657 7432 4 2588 
3 — 4 8004 2652 9428 ~- 5416 
4 5645 2653 5073 229 
5 _ 154 2653 4919 75 
6 — 102 2653 4817 _- 27 





The numerical values of the polynomials in x occurring in the above formulae 
are given in Table ITI. 


IV. EXAMPLES OF THE TYPES OF PROBLEM TO WHICH THE 
Expansion Has Bren APPLIED 


(a) The asymptotic approach to Behrens’ integral 


Fisher (1926) developed the ordinate and integral of Student’s distribution 
in a series of powers of n~’, giving the polynomial coefficients so far as the 
fifth adjustment. The purpose of this expansion was to supply sufficiently 
accurate values of the probabilities corresponding to any values of ¢ for values 
of n beyond the range which it was proposed to tabulate. 

With Behrens’ extension of Student’s test there were even stronger reasons 
for using a similar method. The direct calculations carried out by Sukhatme 
(1938) are very much more laborious than those needed for Student’s integral. 
At any single level of significance, various values are required for three param- 
eters, provided by the two numbers of degrees of freedom of the two samples, 
and the estimated ratio of. the variances of the two means. For functions of 
many variables, there is a great advantage in the use of explicit formulae in 
which the several variables may be substituted, and there is much to be gained 
by extending the use of such formulae over regions too extensive for complete 
tabulation. Finally, it should be noted that the logical situation in which we 
would prefer to rely on the separate estimates of variances from the two samples 
rather than on any process of pooling these estimates, is of more frequent occur- 
rence with large samples than with small, and is particularly applicable to cases, 
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such as arise in Physics and Astronomy, in which we wish to compare estimates 
of the value of the same quantity (a) from relatively ample data of low intrinsic 
accuracy, and (b) from a small series of observations of relatively high precision. 
When, as often happens, the estimates of precision of the means obtained in 
these two ways are of the same order of magnitude, the only satisfactory test 
is that based on Behrens’ solution. The asymptotic expansion is particularly 
suitable for evaluating the percentiles for this special application. There were 
thus four manifest advantages of the asymptotic approach to Behrens’ integral 
(Fisher 1941): 


(i) a check on Sukhatme’s values, obtained by a completely independent 
method, and applicable at least for the higher values of n, and nz , 

(ii) greater accuracy than could be obtained for percentiles from Sukhatme’s 
table for values of n, and n, greater than 12, 

(iii) a wider range of levels of significance in the region to which the asymptotic 
expansion is applicable, 

(iv) the theoretical guidance offered by the algebraic form of the leading 
terms of the expansion. 


(b) The fiducial distribution of the binomial parameter, p 


When discussing the application of the fiducial argument to discontinuous 
observations, Fisher (1959) found that the mean of the fiducial distribution 
of p, the parameter of the binomial distribution, for given observational fre- 
quencies a, b out of N, was 


s.t.2%-* .2~s , Boe... 
Pp — N + ON? 2N® + 2N* (4) 


if x were taken to be normally distributed. | 


On the other hand, the mean of the Bayesian distribution a posteriori, using 
the Bayesian probability a priori 


1 
dp, 
TV pq 


a+} a, b-o@ b-a,;b-a_. 
N+1 Nt on ~ on + oN (5) 


Asymptotic agreement between these means appears when allowance is made 
for the effects of departure from normality in the binomial distribution, which 
are appreciable in the expression (4). 

Direct application of the asymptotic expansion, using the six adjustment 
terms gave the following expansion for the binomial variate in terms of the 
normal deviate x 


a = pN + 2VNopg + dq — pia’ — 1) 


i 3 3 
x } > 2u° — 14% 4 oo? 
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PERCENTILE POINTS OF DISTRIBUTIONS 
which, after inversion, gave* 


ab (b — a)(2x” + 1) 
6N * 
+ (naw + 26ab)2* + (—7N? + Mate tea 
72N*? /ab 


whence, substituting its average value for each power of z, the mean of the 
fiducial distribution is 


b— b— 
“an + oye @ 


agreeing so far as the fourth term with (5). 

The expansion (6) also provides a ready means for comparing the two distri- 
butions with respect to other properties. For example, although the means 
are in agreement, the asymptotic fiducial distribution has the higher variance. 

See also the alternative treatment in Fisher (1957). 


(c) Quantitative inheritance 


Panse (1940) has described a statistical technique for the study of quantitative 
inheritance, in which genetic models, based on data from the F, and F; gene- 
rations, are used to represent the constitution of particular characters. The 
statistical consequences in the population, corresponding to these models were 
assessed, using the cumulant function of the joint distribution of the F, pheno- 
typic value, the mean of the F; progeny, and the genotypic variance of F; 
progeny. These functions provided the data for expressing an attribute of the 
F; progeny in terms of the F, phenotypic values, and thus the effects of selection 
in the F, phenotype on the mean value of the F; progeny could be determined 
by integration over the F, distribution. When the intensity of selection was 
assigned, the limits of integration were calculable from the asymptotic expansion 
of the deviate. 


V. A Cuass or DIstTRIBUTIONS WITH A FINITE CONDENSATION AT ZERO 


The Poisson Series, a discontinuous distribution of positive integers, is well 
known to have the simple series of cumulants 


Ky = ™, 
for all values of r. Correspondingly, the cumulative function is 
K = me‘ — I), 
and the characteristic function 
M = exp {m(e‘ — 1)}. 


* For the remaining terms see Fisher (1959) 
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It is less well known that if, for all values of r, 
x, = rlm, 
or, if a scaling factor be introduced, 
x, = rla’m, 
the distribution, derived from 


mait 
K= 1 — ait 


M sie emi et) 


ot NE ey ene 
=e eval! ati) 


is continuous over the range of positive values, with a finite condensation at 
zero. For (1 — ait)~* is the characteristic function of the Eulerian distribution 


maa 


2 — 
>, ss ’ 


for 2n degrees of freedom, or of the sum of n random variables, each distributed as 


, x > 0, a positive 


ea dz 
a 


Hence the distribution is that of the sum of a number of such variables, when 
the number is distributed in a Poisson Series of parameter m. The variate 
is, therefore, never negative, but is zero with finite frequency 


é 


Over the range of positive values, the distribution is continuous, and can be 


expressed as 
[m. emer (a pe) dx 
va a 


where J, is a Bessel function, specified by 
ites 
2’ .4°.6 


The distribution was first recognized (Bennett 1952, Fisher 1954) as charac- 
teristic of that of the length of germ plasm still heterogenic at an advanced 
stage of inbreeding, but its intrinsic incorporation of a finite condensation at 


Lh) = tut gg + ue +e 
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zero makes it appropriate to a number of natural phenomena, a good illustration 
being the rainfall of an arid region. A continuous model for rainfall is unsatis- 
factory for such regions. For many localities it is preferable to use a model 
ascribing the total rainfall for a given period, for example, a month or a year, 
to a number of showers, the number being a random sample from a Poisson 
series with parameter m, the rainfall of the showers having positive values 
only, representable by the Eulerian distribution 


a dx 
p! 


where p can be small, as in the previous example where it is actually 0, but 
must be > —1. 

The advantage of this type of distribution, relevant for the purposes of the 
water engineer, is that there is a finite probability, namely e~”, of no rain, 
whereas a continuous distribution would make this probability zero, contrary 
to experience. 

For the asymptotic expansion, when m is sufficiently large, we may take 
the exact values for the mean (m) and variance (2m) and for the measures of 
non-normality 


, - 15v2 
mV m 


These yield the six adjustments, to the normal deviate x 


47° — 2 
IV — “384m 
Vv 32‘ + 22? — 11 
480m? -V/ 2m 


2-1 _ 96a" + 164° — 7672 


mV am 46080m° 


The coefficients of Table IV give a rather comprehensive tabulation of the 
distribution, when m is sufficiently large for the accuracy required, and the 
levels of significance of interest. These large values of m would be troublesome 
to use in direct evaluation. For sufficiently small values of m, however, the 
probability that the variate x exceeds any limit X may be evaluated as 

eer eX om 
P =e"e* = — 
pit 2 a! isi iy! 
which may be recognized as the probability that a random Poisson variate 
with parameter m shall exceed a random Poisson variate with parameter X. 
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TECHNOMETRICS 


An Approximation to the Negative Moments of 
the Positive Binomial Useful in Life Testing 


W. MENDENBALL' Anp E. H. Lenman, Jr. 
North Carolina State College 
Raleigh, North Carolina 


The purpose of this paper is to obtain the mean and variance of the maximum 
likelihood estimator of the scale parameter of a Weibull distribution where the sample 
is censored at a fixed time. It will be shown that these moments are functions of the 
negative moments of the positive binomial distribution. A simple approximation 
is obtained for the negative moments of the positive binomial, thus giving an approxi- 
mate expression for the mean and variance of the estimator. 


1. INTRODUCTION 
A family of distributions useful in fitting failure populations is the function 


fi) = () terre where t>0,a,m> 0. (1.1) 


It is known by engineers as the Weibull function as a result of its application 
by W. Weibull to several problems in engineering (10). The parameter m is 
called the shape parameter, the choice of which provides a wide range of distri- 
butional forms including the negative exponential (m = 1). The properties 
of the distribution as well as several methods of estimating the parameters, 
m and a, are discussed by Kao (7). In the case where the shape parameter is 
assumed known and we desire to estimate a, the maximum likelihood estima- 
tion of a reduces to the problem of estimating the parameters of the negative 
exponential distribution due to the invariance property of maximum likelihood 
estimators. Problems of estimation and tests concerning the parameter of a 
negative exponential distribution have been investigated by Epstein and Sobel 
(2, 3) and others. These results are applicable to the Weibull distribution when 
m is assumed known. An extensive list of references is given in (8). 

Data obtained from a life test have the property that the observations are 
received in order of magnitude. The first item failing will be the smallest, the 
second will be the next to the smallest, and so on. We can take advantage of 
this property by concluding the test after either a fixed number of items fail 
or after a fixed length of time has elapsed. Such sampling procedures are called 
censored sampling. By increasing sample size, the length of time for a fixed 
number of units to fail will be reduced. One is therefore able to buy time at 
the expense of increased sample size. Censoring after a fixed number of items 


1Now at Department of Mathematics, Bucknell University. 
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has failed has the advantage of providing a more or less uniform amount of 
information in repeated sampling with the disadvantage that the length of 
testing time varies from test to test. Censoring after a fixed length of testing 
time provides a constant length of test time in repeated tests with the amount 
of information varying from test to test. One advantage of the latter censoring 
procedure is that it simplifies the problem of test scheduling in a production 
process where information from periodic production lots must be obtained at 
regular time intervals. 

In this paper we are interested in investigating the properties of the maximum 
likelihood estimator of the scale parameter, a, where it is assumed that m is 
known and that the sample is censored at a fixed time 7. It can be shown that 
the maximum likelihood estimator of a is 


Ym +m@—n™ 
4 = SO 
r 


number of items placed.on test 
length of test time 
number failing before time T 
number surviving to time T 
; = length of life of the 7th ordered failure (where ¢, is the smallest, etc.) 


In an experimental situation, the experimenter would like to know how 
large the sample size should be in order that the estimates fall within a given 
distance of the true value a. One method of expressing the “distance” numerically 


would be to use the mean square error, E{(@ — a)}, where E represents the 
operation of taking the average or expected value of the square of the difference 
(@ — a). The problem may then be restated, ‘(How large should the sample 
size be in order that E{(@ — a)’} = K, where K is some desired value?” Or, 
we may ask the converse, “For a given sample size and censoring time, what 
will be the size of the mean square error?” 

It will be noted that the mean square error is a function of the mean and 
variance of @&: 

E{(@ — a)*} = V(@) + {E@ — a}? (1.3) 

where 


V(@) = variance of & and 
E(@) = average or expected value of @&. 


We will derive an expression for the mean and variance of @, showing that 
both are functions of the negative moments of the positive binomial, i.e., func- 
tions of E{1/r*} where r is distributed as the truncated binomial, 


oe teen. (1.4) 
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A simple approximation to these negative moments will be derived and com- 
pared with some tabulated exact values for k = 1, 2. 


2. Tue MEAN AND VARIANCE OF @ 


Let ¢ equal the length of time to failure for an experimental unit and let 
the probability density function for ¢ (failure distribution) be 


f(t) = (Bete, t>0 (2.1) 


A random sample of n units is tested until a predetermined time 7’ at which 
time the test is concluded. Times to failure for r observations are observed 
where it will be noted that r is random. Thus we observe ¢, , /. , --- ¢, where 
t; is the jth ordered observation. A total of r units failed and (n — r) units 
survived. 


For the sake of convenience and without loss of generality, we can measure 
time in units of size 7’. Then 


f(z) = ' ee,” were, 


a 


z>0O and B = Fm 


It can be shown that the maximum likelihood estimator of 8 is 


8 = (| Sar +@—n. (2.3) 


t=1 


It would follow that @ = 78. 

In finding the mean and variance of 8, and hence 4, we will restrict our results 
to the set of estimates for which r > 0. When r = 0, the maximum likelihood 
estimate is infinite. In order to obtain more information on a, either n or T 
(or both) should be increased. 

The expected value of @ for the case m = 1 is given by Bartholomew (1). 
It can be shown that the expected value of 8 remains the same, regardless 
of the value of m and that 


ae il 1) . 
HG) = 6-4 + nB(t 1 (2.4) 


p=1-—q and r=1,2,--: ,n. 
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The variance of # is equal to 


via = v(t x 27) +n () + Sn Cor. (2 p> 2", 1) 


(2.5) 


Denote the expectation or variance of a quantity, given r, as E, or V, , respec- 
tively, and let P(j) be the probability that r equals 7. Then considering the 


first term of (2.5), 


vt Sar) =) ve, <1, 
1 


: 
v(2 = e?) = P(1)V(2") + P(g) V2 
oar. 2 


V(x") 


n 


+o + Pe) 
= Ver) D2 Pe) 


1 = 
= a(n 
It can be shown that the variance of x”, given x < 1, is equal to 


Via") = 8 — q/p” 
from which it follows that 


v(t Sat) = @ - ap BU/) 


1 


The third term of (2.5) involves the covariance of 1/r >> 2” and 1/r. 
1 


on (tat t) a a) — a ape 


1 


Note that z(t > 27) = E(x”) 
1 


and hence that x(2 . 2") is 2 Poe.(t p> 2?) = F(z") 


Also, 2 Dat) = Pager) +72 x2 > 27) 


ro 4M Ba) 


Utilizing (2.11), 2% : 21) = B2")6(2) 


(2.6) 
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Substituting (2.12) and (2.14) into (2.10), we find that the 
1g. 1\_ 
Cov (t doa? 5) = 0 
vi) ~ ot Ee) +0) 
1 
- »(!)rery + #0(!) 
= 24) re )+nv S 
-1() +8 — wr 
r UPN; 
Summarizing; the mean, bias, and variance of 8 are 


E(@) = B — a/p + nz) sid 
Bias of 8 = E(8) — B 


“ nB(+) — oe (2.17) 


v() = w'v(4) + @ — a/rx(?) (2.18) 


It follows that 


+ 


E(@) = T"E(8), (2.19 
Via) = T*"Vié), (2.20 
and Bias @ = 7T™ (Bias of A) (2.21) 
The mean square error of # is given by the equation 
E( — 6)’ = V(8) + Bias 8)” (2.22) 
and it would follow from (2.20) and (2.21) that 
E(@ — a)” = T’"E(6 — 8)’. (2.23) 


It is apparent that the mean square errors of § and @ are functions of the 


moments of 1/r since 
rt) - (3) -[)] 220 


and can be computed given n, 7’, an approximate value of a, and the first two 
moments of 1/r. Desk computation of E(1/r) and E(1/r’) is exceedingly laborious. 
To overcome this difficulty, we present an abbreviated set of tables in Section 3 
and a simple approximation formula in Section 4. 
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3. NEGATIVE MoMENTs OF THE PosItTIvVE BINoMIAL DiIstRIBUTION 


Examples of the use of the negative moments of the positive binomial distri- 
bution are given in papers by Stephan (9) and Grab and Savage (5). The latter 
paper discusses in particular the situation arising in sample survey theory 
where it is desired to find the variance of a sample mean, 


co. where the Yi; 7=1,2,---71, 


are identically distributed with mean equal to u, variance equal to o” and r 
is a random variable distributed as a positive binomial variate. The variance 
of 7 is given by 


V@ = B(+), thus requiring 


the first negative moment of r. Reference is made to extensive discussions 
of this problem in (4) and (6). A table of E(1/r) is given in (5) for 


p= .Ol, .05(.05).95, .99 and n = 2(1)20. 
TABLE 1. 
Negative Moments of Positive Binomial Distribution E(1/r*) 
k=l 


15 





ions 
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TaBie 1—Continued 


k=2 


95 


and for 
p = .01, .05(.05).50, .99 and n = 21(1)30. 
In this paper we present tables of E(1/r*) for 
p = .05 (.05) .95,n = 5, 10, 15, 20, 30,40, and 
k = 1,2,3,4. 


Table 1 presents the negative moments for k = 1, 2; Table 2 for k = 3, 4. 
The tabulated results are obtained from the relation 


7 1 ae 1 — z n ro on-r 
»(5) ~i1-¢ La (")o'a 
and are rounded to five decimal places. 


4. AN APPROXIMATION TO THE NEGATIVE MOMENTS OF THE 
Positive BInoMIAL DIstRIBUTION 


Several approximations have been given for negative moments of the positive 
binomial distribution. Stephan (9) achieves an approximation by expanding 
E(1/r) in a series of inverse factorials which converge to the true value of E(1/r) 
as more and more terms are added. An extension of this result is utilized to 
obtain the higher negative moments. Stephan’s procedure permits the calcu- 
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lation of E(1/r*) to any degree of accuracy and enables one to calculate upper 
and lower bounds for the exact value. Its disadvantage is that the formula 
for E(1/r*) lacks computational simplicity. 

Numerous references to approximations for E(1/r) are given in (5). Specifi- 
cally, Grab and Savage mention the approximation 


a(t) wo (4.1) 


and state that one can expect two place accuracy for values of np > 10. We 
will derive an approximation which is nearly as simple as (4.1), which can 
be used to approximate the negative moments for k > 1, and which appears 
to be more accurate than (4.1) for values of np > 5. 

As previously mentioned, Stephan derives his approximation by expanding 
1/r into a series and then finding the expected value of the first ¢ terms. In 
other words, he obtains the expected value of an approximation to-1/r. Rather 
than approximate 1/r, we will approximate the probability distribution func- 
tion of r and then find the exact value of E(1/r) for the approximating distri- 
bution. In consideration of the moments desired, the Beta function is an obvious 
choice. 

The Beta function is given by the equation 


a 21 = . 
f@) > Bla, b) 


(a — 1)!(b — 1)! 
Be, = DT 


Let y = r/n where r is the positive binomial variate. Then the probability 
distribution of y will be a histogram ranging from 1/2n to (2n + 1)/2n. We 
will approximate this function with the Beta function, f(z), although obviously 
an error will be introduced because of the discrepancy between the ranges 
of y and z. This is negligible for large values of np. 

The parameters of the approximating Beta function were obtained by equating 
the first two moments of the two distributions. The two equations are: 


-2,56>0, on8d O0<¢35 1. 


Gah = /- 9) (4.3) 


EEG EOE D 7 mat eV 4) 


where the moments of P(r) are obtained from (9). In order to simplify the 
solution for a and b, let us assume that np is large and hence that (1 — 9°) 
is approximately equal to one. The values of a and 6 are then the simultaneous 
solution of 
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na__», 
— 


n'a(a + 1) x 22 
atiatitn tt" 


a=(n—1)p 
b=(n— lg 


The approximating distribution is then 
jo) = a - a 
Bia, b) 


(n—-1)p-1 (n—1)qa-1 
_¥ (1 — y) 
Ba, b) where 0< y <1. 


The negative moments of y are easily obtained from (4.5) yielding 


a) Mic 


_@-—k- Ib - @ kb — Wt [a — Diep — 1)b — 1)! 
(a+b—k-—1)! (a+b— 1)! 


_@+b—-NG+b—-2---G@+b-B 
(a — 1)(a — 2) --- (a—k) 


Specifically, from (4.7), 
. We y c 


ously 
anges 


1\_ @—2)™@ — 3) 
(4) “@-te- 9 (4.9) 


lating 
As noted in (4.1), in order for the moments (4.6) to exist, it is necessary that 
(4.3) a=(n—Il1)p>k. (4.10) 


The approximations to the moments of 1/r are easily obtained from the 
moments of 1/y. Since 1/r = 1/(ny) 


25) iw + (4) (4.11) 


(2) S a S 5 


In particular, 
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> 


(1) _ _(n — 2)(m — 3) 
x ) ™ n(a — 1)(a — 2) 


r 
Utilizing (4.12) and (4.13), it can be shown that 


1\ __ @ — 2)(n — a — 1) 
v(?) <"n@ — 1" — 2) to 


In order to distinguish these approximations from others, we will call them 
the Beta approximations. 

The accuracy of these approximations will be discussed in Section 5. It might 
be well, however, to comment upon the errors introduced by our simplifying 
assumptions. The first error was caused by the error in approximating the 
range of the histogram of P(r) as given by equation (1.4). The second was 
caused by the assumption that gq” is small and hence that the equations for 
the solution of a and b are approximately as given in (4.3). Both assumptions 
were made to obtain simplicity in the approximation and both errors diminish 
as np increases. We will see in Section 5 that the approximation is quite good 
for moderate sized values of np in spite of the assumptions. 


5. A CoMPARISON OF THE APPROXIMATION WITH TABULATED 
Vauues or E(1/r) anp E(1/r’) 


A comparison of the approximation to E(1/r) with the tabulated values 
for various combinations of n and p is presented in Table 3. In addition, the 


approximation (4.1) is included. Proceeding from top to bottom, the three 
values in each cell are the tabulated value, the Beta approximation (4.12), 
and the approximation (4.1), respectively. All values are rounded in the sixth 
decimal place. The tabulated value for n = 100, p = .90 is omitted due to 
computing difficulties. Also, no approximation is given for the cells where a < k. 

Examining Table 3, we note that the Beta approximation is correct to two 
decimal places for np > 5 and in several instances is correct to six decimal 
places. The approximation (4.1) also appears to be quite good although the 
Beta approximation would, in general, appear to be more accurate and hence 
more suitable for use in computing the mean square error of the estimator 
& mentioned in Section 2. As n and p increase, (a — 1) — (n — 2), g — 0, and 
as a result, both the Beta approximation and the approximation given by 
equation (4.1) approach 1/n. 

Stephan (8) gives two examples of the use of his approximation. In the first 
example, n = 100, p = .1 and E(1/r) = .111527, correct to six decimal places. 
Utilizing two terms of the series, he obtains E(1/r) ~ .108675; for three terms, 
E(1/r) = .110548. The Beta approximation is .110112. In the second example, 
n = 1000, p = .3 and E(1/r), correct to nine places, is .00334116. Two terms 
of the series yield E(1/r) = .00334108 and for three terms, E(1/r) ~ .00334116. 
The Beta approximation is .00334114. Thus, in these two examples, it requires 
three terms of Stephan’s series to equal the accuracy of the Beta approximation. 

The tabulated value of E(1/r’) and the Beta approximation are given in 





alues 
, the 
three 
1.12), 
sixth 
ue to 
ee fe 
0 two 
cimal 
h the 
hence 
mator 
0, and 
en by 


e first 
places. 
terms, 
ample, 
terms 
34116. 
aquires 
nation. 
ven in 


AN APPROXIMATION TO NEGATIVE MOMENTS 


TABLE 2 
Negative Moments of Positive Binomial Distribution E(1/r*) 


k=3 
15 
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TaBLe 3. A Comparison of the Tabulated Value of E(1/r), the Beta Approximation, and 
the Approximation Given by Equation (4.1) 
Entries: Tabulated Value 
Beta Approximation 
Approximation (4.1) 


20 30 


.576912 . 428114 .323807 . 111527 

1.00000 .491228 .327586 . 110112 

. 909091 . 476190 .322581 . 109890 

.416735 . 194724 . 122147 .089094 -034159 
-470588 . 191489 . 121212 -088785 .034146 
. 434783 . 188679 . 120482 088496 -034130 


- 229109 . 105990 069152 ‘ 051356 .020206 
. 228571 - 105882 -069136 .051351 - 020206 
. 222222 . 105263 . 068966 .051282 .020202 


. 150568 .073123 .048345 .036116 -014348 
. 150943 .073171 048359 .036122 .014348 
. 149254 .072993 -048309 .036101 .014347 


. 112523 -055885 .037180 .027857 
. 112676 .055901 -037185 -027859 .011124 
. 112360 .055866 .037175 027855 .011123 


TaBLeE 4, A Comparison of the Tabulated Value of E(1/r3) and the Beta Approximation 
Entries: Tabulated Value 
Beta Approximation 


20 30 100 


- 252093 . 149937 .014559 
.049123 . 159483 .013520 


.017077 -008627 -001198 
.016282 -008467 -001196 


.004985 -002716 -000413 
.004978 .002714 .000413 


. 125262 4 : 002375 .001320 . 000207 
. 166667 ‘ 7 -002378 -001321 .000207 


054698 . ‘ -001388 -000778 
.057692 ‘ = .001389 .000779 .000124 





114559 
13520 


101198 
101196 


100413 
004.13 


00207 
100207 


00124 


——_— 
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Table 4, the tabulated value followed by the approximation. The agreement 
appears to be quite good for moderate sized values of np. 


6. AN EXAMPLE OF THE DETERMINATION OF THE MEAN SQUARE ERROR OF @ 


We plan to conduct a life test in which time, 7’, and sample size, n, are dictated 
by the amount and availability of material and equipment. We desire a measure 
of the goodness of estimation which the money represented by n and 7 will 
buy. As mentioned in Section 1, the mean square error of the estimator provides 
this measure. 

For example, let us assume that we wish to estimate the parameter a of 
the Weibull distribution for the estimation problem discussed in Sections 1 
and 2. Further, let us assume that from previous experience, we know that 
the failure distribution has a form represented by the Weibull distribution 
with shape parameter m = 2 and that the mean life lies somewhere in the range 
of 80 to 120 hours. We might choose 7 = 200 hours and sample size n = 40. 

In calculating the mean square error of &, we will use equations (2.17, 2.18, 
and 2.22) to obtain the mean square error of 8. Utilizing (2.23), we can deter- 
mine the mean square error of @. 

Looking at equations (2.17 and 2.18), we see that we need to know the values 
E(1/r), V(1/r), 8, p and q. Since 6 (or a) is unknown and is the parameter 
we wish to estimate, the best that we can do is to substitute a guessed value 
for B. 

The mean value of ¢ for the Weibull distribution is 


E(t) = aver(1 + 1) (6.1) 


Pe (6.2) 


+) 
m 
Since we expect the mean life to fall in the range 80 to 120 hours, we might 


choose E(t) = 100 as a guessed value. Then, for our example, the guessed 
value of a would be 


100 | 
1 
t= 
(+1) 
Tables of I(x) can be found in most standard handbooks of mathematical 
tables. We note I'(1.5) = .88623. Then 


~ (_100_\ 
ams (00) —? (6.3) 


From (2.2) and (2.4) 


(6.4) 
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q = 6 '/F = @7 17-8180 _ 194390 (6.5) 


and 
p =1-—q= .95680 (6.6) 


The values of E(1/r) and V(1/r) can be determined from the Beta approxi- 
mations (4.12) and (4.14). 


(4) = oo . where a = (n — 1)p = (89)(.9568) 


r 
= 37.315 


fio nap aT 
x(?) = 40(86.315) ~ -°262 


-{1\ __ @ — 2)(m — a — 1) 
v(2) ~ na — 1)?(a — 2) 


_ (88)(40 — 37.315 — 1) 
~ - (40)*(86.315)*(35.315) 


Substituting the values into (2.17) and (2.18), Bias 8 = nE(1/r) — g/p — 1 


0432 _ 
9568 


= .00000085927 


= 40(.0262) — 1 


= .0013 


vi) = n'v(2) + @ — a/pe(4) 


.0432 


= (40)*(.0000009) + (.318)’ — (9568)? 


(.0262) . 


= .00284 


The mean square error of # and @ can be calculated from (2.22) and (2.23) 
respectively. 


E{(6 — 6)"} = V(8) + (Bias 8)’ 
.00284. + (.0013)? = .002849 


In this case, the square of the bias is so small compared with the variance 
that it can be neglected. 
E{(@ — a)"} = T’"E{(6 — 6)"} 
= (200)*(.00285) = 4,574,000 
When the bias equals zero, the mean square error equals the variance and 
consequently, for this case (bias approximately equal to zero), the standard 


deviation of the estimator, &, is 2139. 
The mean square error as a function of 8 is presented in Figure 1 for n = 10, 
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20, 40, and 70. For the above example, we read from Figure 1 (n = 40) the 


mean square error of 6 is approximately .0030 for 8 = .318. 


The authors wish to thank Dr. A. H. E. Grandage for assistance in the prepa- 
ration of Tables 1 and 2. 
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Order Statistics from the Gamma Distribution 


SHANTI S. GuPTA 
Bell Telephone Laboratories, Incorporated 
Allentown, Pennsylvania 


This paper deals with order statistics from a gamma or x? (Pearson Type III) 
distribution. Expressions are derived for the moments of an order statistic and for 
the covariance between two order statistics. A table of moments (about the origin 
as well as the mean) is presented. Equations are derived which are used to solve for 
the percentage points of the order statistics and a table of the percentage points is 
given. The problem of the best linear unbiased estimate of the scale parameter of 
the distribution which is based on order statistics is, briefly, discussed for the cen- 
sored (type II) case. The modal value of an order statistic is derived and a table of 
these values is given for the cases of the smallest and the largest order statistics. 
Applications to life-tests, extreme values, reliability and maintenance are described 
and illustrated in some cases. 


1. INTRODUCTION 


Suppose that we have N independent random variables x, , z, --+ , Zy each 

distributed according to the gamma probability density function (p.d_f.) 
—_ 

(1.1) g(x) =e TO’ 0<r<o 
where r is assumed to be a positive integer. The distribution (1.1) arises, for 
instance, as the convolution of r independent and identical chance variables 
with the exponential density function e~*. Let G,(x) be the cumulative distri- 
bution function of the z, then it is easy to see that G,(x) can be written as a 
partial sum of probabilities in a Poisson distribution as follows 


(1.2) G@=> os. 


Let x) denote the kth smallest observation when the z’s are arranged in the 
increasing order, viz. 


(1.3) Lay S 2) ¥e's Lay Sess SX), 


and let y = rq) , then we are interested in the moments and distribution of y. 


2. MoMENTs OF THE kth OrpErR STatTIstTIc 


The ith moment about the origin, u{(k, N) of the kth order statistic from 
a sample of size N from the gamma p.df. g,(x) is given by 
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, T be aac cect ee 
uk, N= Expl w bp! 


© ;  r-1 “wy oa 
f E hear a. Te 


which reduces to (details omitted) 


, sas I ek So fb — 1 
wk, N) = G—prw — ir 2 ' m4 P ) 


(2.1) 


p=0 


(2.2) 
(r-1) (N-k+p) T(r +it+ m) 

m ev ee k Fee ee NT TOT 
ao +) W—k+p +) 

r—1l gi 

where a,,(r, p) is the coefficient of ?” in the expansion of (Ss 


i=0 J: 


"For k = l, 


we have from (2.2) 


awe a rr +i +m) 
(2.3) mi, N) = Foy De aml, N — 1) ae 


The zth moment of the kth order statistic can be expressed in terms of the 
ith moments of Ist (smallest) order statistic as follows 


; i N! Ss _elk — 1)\ ni, N-k+p+4+)) 
24) Wik, §) = Era & ( *( p ) (N-E+p+h) 


Also, the coefficients a,,(r, p) satisfy the following basic recursion relation. 


a,(7, Pp) = an(r,p — 1) + Gnalr, p — 1) 
(2.5) 1 1 
+ 5j m2", BP — Wee + @— Dl amr P = 4) 


Since for p = 1, these coefficients are 1, 1, 1/2!, -++,1/(r — 1)!, we can compute 
them for any value of p by repeated applications of (2.5). 
If we define 


Hk, N) = [GGL — G,@)"*g,() dz 

(2.6) . 
_ & = 7 As k)! T(r) uli, N), 

then v{(k, N) satisfy the following recursion formula 

2.7 veh, N) = (—9"(% > ye + ak +0). 

A slightly more generalized form of (2.7) is the following one. 

2.8) rie, N) = (12k + pV — 8 +2) 


where £ is any positive integer < N — k. Also, the relation (2.4) can be written as 


(2.9) vi(k, N) = > —" 3 Ned N-—k+p+). 





ORDER STATISTICS FROM THE GAMMA DISTRIBUTION 
From (2.8), we have, by putting 8 = 
(2.10) i(k, N — 1) = vf(k, N) +k +1,N), kSN—1. 


Table I gives the first four moments (about the origin as well as about the 
mean) of the kth order statistic for selected values of k, N, and r. 


TABLE I 


The entries in this nage are the moments about the origin yy (k,N) 
of the eh order statistic of a sample of 
size N from a standardized gamma distribution with parameter r. 


pel 


2.00000 | 00000 2.00000 | 00000 6. 6.00000 | | 24.0000 0000 i. 1.00000 —*| 2.00000 CY 00000 9.00000 CO 00000 


wnr 


Fune 


UI we 


\ 


3 
4 


wr 


UI UNH 


BHR HH OWOMIANLSWUNH WMd 


Be Ree 
UU e 


«500000 


1.50000 


. 233333 


200000 
450000 
783333 
28333 
28333 


-0769231 
-0714286 
.0666667 


«500000 
-50000 


- 222222 
-05556 
- 72222 


- 125000 


-0200000 
0669136 
-150941 
.287789 
-503002 
-841253 
-38908 
34171 
-27069 
.1286 


.0165289 
-0138889 
.0118343 
-0102041 
-00888889 


«750000 
«2500 


222222 
80556 


-9722 


+0937500 
«607639 


+00600000 
-0283045 
-0549075 


-208246 


00347222 
:00273100 
00218659 


-00177778 


1.50000 
46.5000 


67.7963 


-0937500 
+ 903935 

6.910838 

88:0914 


-0384000 
- 315150 
1.78711 
10.3267 
107.533 


0185185 
°137807 
- 669835 
2.90439 
14.0379 
126.232 


-00999583 


0696546 
-308189 
1.15203 
4.21866 
17.9656 
144.276 


00585938 


-0359511 
.161765 
-552229 
1.75183 
5.69876 
22.0545 
161.736 


-00365798 


0234705 
-0931325 
- 299030 
- 808728 
2.45831 
7.31898 
26.2647 
178.670 


-00240000 | 


-0149798 
-0574336 
opfoaae 
«482929 
1.25453 
2-82 
9.05819 
30.5663 
195.126 


+00163923 
-00115741 


-000840307 
-000624740 
-000474074 


+ 250000 


1.25000 


-111111 
. 361111 
1.36111 


-0625000 

-173011 

«423611 
1.42361 


-0400000 
. 102500 
213611 


1.49139 


.0204082 
0481 859 
881859 
*150686 
-261798 
«511796 


-0156250 
0360332 
206 38109 


:277422 
» 527422 


235i51 


«0100000 
pan 
0379707 
0583790 
0881352 


-00826446 
00694444 
-00591716 
-00510204 
00444444 


+ 250000 
2.25000 


-O740741 
«324074 
2.32407 


-0312500. 

+ 105325 

355324 
2.35532 


-0160000 
-0472500 
-121324 
«371324 
2.37132 


.00925926 

10252593 

-0565093 

«130583 

- 380583 
2.38058 


-00583090 
0150902 
0310902 
.0623401 
136414 
° 386416 
2.38641 


.00390625 
00972726 
0189964 
0349965 
30062462 
-140321 
1390323 
2.39031 


.00274348 
-00664373 
-0124806 

0217399 
-0377396 
0089337 
«143057 
«393073 
2.39306 


«00200000 
°00474348 
~00804973 
0144806 


2.39504 


-00150263 
200115741 
-000910332 
.000728863 
-000592593 


-562500 
11,0625 


111111 
. 840278 
1.20069 


0251383 


13.9108 


-0144000, 
0645563 
-244001 
1.12692 
12.9086 


- 00694444 
*028 0111 
-0885840 

286547 

1.21133 


13.1595 


-00374844 
-0140942 
.0400589 
108285 
319854 
1.27504 
13.3458 


.00219727 

.00785897 

.0208089 

-0505235 

.124609 

» 346592 
1.32522 
13.4898 


-00137174 
-00472642 
-0118998 
.0269073 
-0595857 
«138302 
«398531 
1.36564 
13.6043. 


-006900000 
-00301248 
-00730466 


1.39885 
13.6977 


-000614712 
-000434028 
-000315115 
-000234277 
-000177778 
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TABLE 1 (Cont. ) 


The entries in this table are the moments about the origin Hy (KN) 
and the central x nts iG) of the «rh order statistic of a sample of 
size N 
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TABLE I 


Cont. 


The entries in this table are the moments about the origin 
ny (k, N) and the central moments oar of the _ order statistic of 
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TABLE I (Cont. ) 


The entries in this table are the moments about the origin uy (k, N) 
and the central moments Hy (K.N) of the k*" order statistic of a sample of 
size N from a standardized gamma distribution with parameter r. 


r=4 
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TABLE I (Cont.) 


The entries in this table are the moments about the origin Hy (KN) 
and the central moments Hy («,N) of the Kk’ order statistic of a sample of 
size N from a standardized gamma distribution with parameter r. 


r=5 
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TABLE IT 


The two entries in each cell of this table are the modal values of the small- 
est (top) and the largest (bottom) of a sample of size N from the standardized gamma 
distribution with parameter r. 
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3. CovARIANCE BETWEEN THE mth AND nth OrpzER Statistics 


Let & = 2; and 7 = 2;,; be the mth and nth (n > m) order statistics from 
the gamma distribution (1.1). Then, we have 


G1) Be) =c f [ ” eG?" 1G.(n) — 6.) 


- [1 — G(n)]"""g-(€)9-(n) dé dy 
where 
ap 
~ (m—1)!(n— m— I!I(N — n)! 


The expression on the right side of (3.1) can be written as 


Be) = 0 D(-y(™ > Tn mB) 


© is r-1 u’ a-1 -— o ha r-1 v' s-1 i 
‘[? [= | rom, ? [= ‘i re” 
whereg = a+n—m-— Bands = N—n+68+1 and the first summation 


on the right side is over the positive integers a, 6 (0 S< a S$ m — 1, 
058 Sn —~™m -— 1). As before, let a,(r, p) denote the coefficient of u‘ in the 


expansion for 
o=5 ui P 
then (3.2) simplifies to 


Fe) = Gey SA" SO’) 


ee a(r,s — 1)T(t+r+ 1) 


t+r+1 
t=0 $s 


(3.2) 


PS (s')\ REY PG 4+r+pt Dal, q— o| 
| (4) X (s+ oo 


(3.4) —" wor a me - pee” z ‘\(" $ ' r ") 


[a $ — Ia lr, g— NT(t+r+)IG+r+p+ 0] 
j! str ittg 4 —" 


From (3.4) we can obtain the covariance (¢, ») by subtracting from E(é) the 
product E()E(n) which we get from the results of Section 2. 
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4. Some REMARKS ABOUT THE DISTRIBUTION OF THE kth OrnpDER STATISTIC 


a. The Mode of the kth Order Statistic 


The modal value ¢ of the kth order statistic from a sample of size N from 
the gamma distribution with parameter r is given by 


(4.1) 5 ewer, Ole — vir, 0) = 0 


where 
r-1 ti 
ag Ss > 
(4.2) V(r, t) py» j! ? r = 1 
= 0, otherwise. 


The function y(r, ¢) satisfies the following difference-differential equation 


(4.3) oun) = ye, ) = ve - 1,0. 


Using equation (4.3) in (4.1) we obtain after some simplification 
(4.4) (—Nt T= 1)y¥(r, 3) a vr, t)] + (N tei k)tv(r =e 1, t)|e’ -* vr, t)] 
+ (k — 1tve, dle’ — yr — 1,9] = 0 


or 


(4.5) ( -—t—- (Ee ~ A = a | — ke’ — (N — 1) 5 


i=0 J: imo J! 


i=0 


Equation (4.5) can be used for an iterative procedure to solve for ¢. The cases 
k = 1 and k = N are of particular interest and we obtain 


(4.6) 


(4.7) 


which give the modal values of the smallest and the largest order statistics, 
respectively. It should be noted that the present forms (4.6) and (4.7) are 
helpful in that the value of the left side expression can be obtained from a 
table of cumulative Poisson distribution. 

It is easy to show that there are unique solutions for (4.6) and (4.7) in the 
open intervals 0 <¢<r—1landr—1<t< ©, respectively, provided N 2 2. 
For N = 1,1 = r — 1, since in this case the mode is the same as that of a gamma 
distribution. Also, ¢ increases monotonically with N for the largest order statistic 
and decreases monotonically with N for the smallest order statistic. 

Another particular case of interest is r = 1, since in this case the gamma 
distribution reduces to the exponential distribution. For this case equation (4.5) 
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reduces to the simple and explicit form 
(4.8) t = log, N — log. (N —k +1). 


Equations (4.6) and (4.7) were solved by iterative methods on IBM 650 
for selected values of N and r. The results of these computations are given 
in Table II. 


b. Cumulative Distribution Function of the Order Statistics 


The cumulative distribution function F,(y; k, N) of the kth order statistic 
y = Za) from the gamma distribution is 


Fak) = Gia pi | SOU - Gea.) de 


(4.9) Fe N! : Hii se Saab 
“ta. Se 
= Tey(k,N —k+ 1), 
where I¢,,,) (k, N — k + 1) is the ratio of the incomplete beta function to the 
complete beta function and is tabulated in [5]. For a given value of y, we can 
compute the c.d.f. F,(y; k, N) by first consulting a table of the cumulative 


distribution of the Poisson or gamma distribution to get G,(y) and then con- 
sulting the table of the cumulative of the beta distribution. Again, 


N! (pt (n - [Sey'\" 
(4.10) Pusk,N) = gaye SG a hE _ 


(4.11) Fy; k, N) = oe b Py; i, a(x i. - He cy | 


The function F,(y; j, j) is the c.d.f. of the largest of jxz’s and using K. Pearson’s 
notation [6] for the incomplete gamma function can be written as 


4.12) Fi, = 61) = (7-1) = {per} 


Also for k = 1 and k = N i.e. for the distribution of the smallest and the largest 
respectively, we get the following simple forms 


(4.18) F,(y;1,N) =1— [1-@,Q@P =1—- “eel 


(4.14) F,(y; N,N) = Gy) = mse ! a 
The c.d.f. of the median i.e. of 2im+1) When N = 2m + 1, is 


(2m + 1)! £ 1)! ( 1)* ("hs eu 
4, = 
( 15) Fy; m+ 1, 2m+ 1) (mI ro | 2d i! . 
Probabilities in (4.10) ... (4.15) can be obtained with the help of a table of 
cumulative Poisson distribution [3] and also with the help of a table of incom- 
plete gamma function [6]. 
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TABLE III A 
r=1 


This table givesthe a-percentage points of smallest order sfatistic 
of size N [1s k <«N,N.= 1(1)10, 25, 50,100, 200, 500, he, BO, and ‘selestea va tugs 


of a) from a standardization gamma distribution with parameter Yr. 


ay 201 05 .10 +25 -50 +75 +90 . -99 


-01005  .05129 «= «.10536 = «28768 = «69315 1.38629 2.302592. 4.60517 


-00503 -02565 -05268 14384 «34657 -69315 1.15129 3 2.30259 
+ 10536 -25310 -38013 -69315 1.22796 2.01014 2.96973 -676 5.29580 


00335 -01710 -03512 P 23105 -46210 + 76753 ‘ 1.53507 
-06071 - 14543 +21791 J 111978 1.63066 : 2.83186 
24263 -45950 -62392 994 ‘ 2.39207 3.36648 - 5.70044 


-00251 senate -02634 ‘ Z es -57565 4 1.15129 
. on 15380 é 4873 -78456 1.13800 1. 1.95992 
"2562 . 38634 $ 95% -41461 =1.94799 ‘ 3.17011 
-64031 - 82631 3 1.83822 -66794 3.64981 $ 5.98771 


.01026 02107. -13863  .27726 = «46051 Six 94209 
.07952 «11905 is -37660 60547 «= 687681 =—«:1.07122 +=—:1.50476 
.20981 .28321 =. -69315 1.02321 1.39983 1. 2.24772 
41945 = .53804 1.15897 1.64114 2.18712 2, 3.42093 
.79091 = .99685 si. 2.04446 2.88396 3.87031 4. 6.21062 


-00855 -01756 4 : +23105 - 36920 7 -76753 
06491 -0971 : ‘ + 4934 -71400 3 1.22312 
.16624 -2242 é P 1.09904 304 1.75406 
-31655 -40525 . 1.21429 1.60490 : 2.46829 
-54163 -67272 . s - 825; 2.37992 - 3.62073 
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For given N, k, (1 s k s N) and a the entries in the above table are the values 
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to the p.d.f. and c.d.f. of the standardized gamma chance variable, respectively. 
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TABLE III B 
r=2 


This table gives the a-percentage points of the «eh smallest order statistic 
of a sample of size N [1 « k < N, N = 1(1)10, 25, 50, 100, 200, 500, 1000, and selected 
values of a] from a standardized gamma distribution with parameter r. 
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For given N, k, (1 ¢ = N) and a, the entries in the above table are the values 
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the p.d.f. and c.d.f. of the standardized gamma chance variable, respectively. 
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TABLE III -C 
r= 3 
This table gives the a-percentage points of the «rh smallest order statistic 


of a sample of size N [1 s k s N, N = 1(1)10, 25, 50, 100, 200, 500, 1000, and selected 
values of a] from a standardization gamma distribution with parameter r. 
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to the p.d.f. and c.d.f. of the standardized gamma chance variable, respectively. 
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c. The Percentage Points of the Order Statistics 


If B,(k, N — k + 1) denotes the a-percentage point of the beta distribution, 
then y,. , the a-percentage point of the kth order statistic is the solution of 


(4.16) G.(ya) = Balk, N —k +1). 


Equation (4.16) can be written as 
r-1 i 

4.17) —va + log. ($9) = tog, 1 — Balk, N — + 0) 
j=o0 9! 


and this equation was used to find y, by an iterative procedure. The values 
of B,(k, N — k + 1) were obtained from [4]; for those values of kand N —k+ 1 
for which B,(k, N — k + 1) was not available, inverse interpolation was per- 
formed in the tables of the incomplete beta function [5]. The latter is not neces- 
ary for k = 1 and k = N since for these cases we have, respectively, 


(4.18) Gy.) = 1 — (1 — a)" 


and 


(4.19) G(y.) = (a)"”" 


5. Best LINEAR UNBIASED ESTIMATES 


We consider the gamma p.d.f. in the form 


— 272/80 zt 
(5.1) g(x, 0) =e gT@ ’ s>¢, @>0 
where @ is to be estimated (r assumed to be known). The best linear unbiased 
estimate of @ based on the ordered observations is well-known for the case of 
a complete sample and is given by 2a/r. However, if we have a censored sample 


(5.2) Lies) BS Voary S ett S Le) > 72 2 (2,7, both known) 


then the problem reduces to determining the vector of constants, b = (b, ,b,,--- 
b,,-r,41) such that 


(5.3) b’ >> b = minimum 


z(k) 


and 
(5.4) b/Ex«) = 6 


where >-;,,, denotes the covariance matrix of the ordered x’s in (5.2) and 
<{,) denotes the row vector (2-,) , Zr,41) *** Ler,))- Since Exu, = 6 wi (k, N), 
the constants b, , b, , «++ , b,,-,,4: are the solution of 


j=f +41, a 
(5.5) 


b,ui(k, N) =|] 





SHANTI S. GUPTA 


s fe-ritl 
j=’ Sb — af 2 baui(k, N) — i 
Z(k) = 
The constants b’s are given by 
b,, 0 
(5.6) Prstl og]? 
b,, 0 
B 1 


where S is the matrix of coefficients on the left side of (5.5). For example, 
ifr, = 1,r, = 2, N = 3, then b, and b, are given by 


? _ wi(1, 3)]~ 
12 2 


_wi2, 3)| {0 
2 


2 


11,8) —_wa(2, 3) 0 


where ¢:: , Coo , Ci2 are the variances and covariance of t1) , 22) in a sample 
of size 3 from a gamma distribution. 

Some constants b’s for the cases of r = 2, 3 will be published later. For r = 1, 
i.e. for the case of the exponential distribution these have been obtained in [8]. 


6. AppLicaTions To Lire-TEsTING, RELIABILITY AND 
OTHER PRoBLEMS AND SoME FuRTHER REMARKS © 


a. Time to rth Failure 


In the following we shall assume the exponential density (5.1) and denote 
by 2,) the time to rth failure in a life-test experiment with replacements and 
by 1 the number of units on test. Then x = nz;,)/@ has the p.d.f. g,(x) of (1.1). 
Thus, the results of this paper will be useful in finding the probability that 
the smallest (or any other ordered value) time to rth failure in N life-test experi- 
ments (with replacements) with n units in each will at least be as large as a 
given value. Also, the expected values and other moments of the ordered times 
to rth failure can be found with the help of Table I. 


Illustration 1 


Life-tests (with replacements) are to be carried on with n = 10 units each 
of the same type, in N = 3 test sets. We would like to know the probability 
that in none of the 3 sets four failures will take place in less than c = } times 
the average life-time. The required probability is P = 1 — F,(2; 1, 3) which 
we can find from a table of the Poisson distribution (the parameter of the Poisson 
distribution being nc = 2) as indicated in formula (4.13). Thus, P = (.857123)* = 
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TABLE III D 
re=4 
This table gives the a-percentage points of the rh smallest order statistic 


of a sample of size N [1 s k s N, N = 1(1)10, 25, 50, 100, 200, 500, 1000, and selected 
values of a] from a standardized gamma distribution with parameter r. 
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For given N, k, (1 ¢ k 4 N) and a the entries in the above table are the values 
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TABLE III E 
r=5 


This table gives the a-percentage points of the rh smallest order statistic 
of a sample of size N [1 s k < N, N = 1(1)10, 25, 50, 100, 200, 500, 1000, and selected 
values of a] from a ‘standardization gamma distribution with parameter r. 
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.6297. Also from Table III, we can put one-sided or two-sided confidence limits 
on the smallest (or any other) time to rth failure in terms of the average life 
or vice versa. For the present example r = 4 and from Table IIID we find 
that the 90% lower confidence limit on x) , the smallest time to the fourth 
failure is .1209 times the average life 6. Again the 90% lower confidence limit 
on 6 based on 2,1) , say, is (10/3.8458) 24) = 2.6002 x, . 

Again by using Table I we find, 

Expected value of the smallest time to the 4th failure = .2449 av. life 

Expected value of the median time to the 4th failure = .3820 av. life 

Expected value of the largest time to the 4th failure = .5731 av. life 

If we have N systems and the rth failure of a unit in any of these means that 

the system becomes unsatisfactory or unreliable (replacements are made till 
the (r — 1)th failure), then the probability that 0, 1, 2, --- N of the systems 
arrive at an unsatisfactory state before a certain time can be obtained from 
the formulae of Section 4 and also individual and joint confidence limits can 
be placed on the reliability of these systems. 


b. Total Life 


If t, , t2, --- , t, are the times to Ist, 2nd, --- , rth failures in n units on test 
(without replacement) that follow the exponential law f(t) = Xe’, then the 
total life 7’ defined by 


(6.1) T=iitéte:-+4ta-—ne, 
has the gamma density function 


Tor 7. 1 


(6.2) g(f) = Ne", 


so that results of this paper will be applicable in comparing two or more popu- 
lations using total life as the criterion. 

The distributions of the statistic T when ¢’s come from exponential and uniform 
distribution are derived in [1] and [2], respectively. 


T>0 


Illustration 2 


The most probable value of the largest order statistic is of importance in 
the study of breakdown voltages of paper capacitors of several layers. It is 
assumed that the size of the conducting particles (impurities) follows an ex- 
ponential distribution and that the particles of one layer are lined up against 
those of the next layer. Under this assumption the total length of the r con- 
ducting particles follows a gamma distribution. We are then interested in the 
mode of the largest of a sample of N from a gamma distribution. 

Let N = 10 be the number of particles per square unit area on each of the 
r = 3 layers. Then from Table II, we find that the modal value of the largest 
length of the conducting particle is equal to 5. 448/ \ where J is the parameter 
of the exponential distribution \e~™’. 
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c. Relation to Order Statistics from x’-Distribution 


The random variable $x3, with 2» degrees of freedom has the gamma density 
function (1.1) with parameters r = v. If y and y’ denote the kth order statistic 
of a sample of size N from the standarized gamma distribution with parameter 
r and x3, (with 2r degrees of freedom) respectively, then we have 


(6.3) Ply’ <¢] = Ply < c| 


(6.4) Ya. = Wa 


where y/ and y, refer to the a-percentage points of y’ and y respectively. Also 


(6.4) E(y’)’ = 2‘Ey’. 


d. Application to Distribution of Ordered Sample Variances from Normal 
Populations 


If s? is the sample variance based on v degrees of freedom from a normal 
distribution with variance o”, then vs?/o” is distributed like a x? with v degrees 
of freedom. This enables us to obtain the moments and distribution of the 
ordered sample variances for y = 2r degrees of freedom from the results and 
tables of the present paper. 
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TECHNOMETRICS 


Parallel Fractional Replicates 


CuTHBERT DANIEL 
New York City 


In studying the effects of many factors on several responses the experimenter is 
sometimes able to say which factors are likely to be influential for each response. 
This knowledge permits more efficient planning since certain effects and interactions 
are now negligible. Some 2?-¢ fractional replicates are constructed that take ad- 
vantage of this knowledge. 


Most industrial products have more than one important property. Experi- 
ments to measure the influence of process variables on several properties are 
often carried on simultaneously. This is likely to be the case when a large part 
of the cost of experimentation is spent in producing material for tests. Thus 
even small batches of a metal, of a plastic, or of a new food product may be 
more expensive to prepare than to test for several properties. As a simple example, 
imagine that three factors affecting the flavor (y,) and texture (y.) of a new food 
product are being studied. An essential oil (A) is to be added that may modify 
the flavor (y,) although it cannot influence the texture (y,). A wheat product (B) 
that may influence the flavor and texture is also varied. Finally, a bran-like 
material (C) without flavor but likely to influence the texture is to be the third 
factor. This state of affairs may be abbreviated as “influence pattern” A-—B-C. 
or as 1-1-1. The first term in such an abbreviation indicates the number of 
factors that are taken to influence only y, . The second term numbers the factors 
that influence both y, and y, . The third term shows the factors that influence 
Y2 only. 

The familiar half-replicate, 2°’, with treatment combinations (1), ab, ac 
and be can be viewed as a full 2” factorial experiment on A and B for studing 
their effects and interaction on y, . At the same time this plan is a full 2’ on 
B and C and so will give information on their effects and interaction on y, . 
Six estimates are obtained then from four “runs’’. If the assumptions made 
about the influence pattern are correct, then effects due to C, AC, and BC 
do not exist for y, , and effects due to A, AB, and AC are zero for y, . There 
is then no confounding in this fractional replicate. 

Proceeding to slightly more realistic examples, assume that the three factors 
A, B and C are known to influence y, , and that three other factors D, E and F 
influence y, . Assume further that each set of factors influences only one re- 
sponse. The influence pattern is now written 3-0-3. The obvious plan would 
be to adjoin a 2° on A, B and C to one on D, E and F. The combined set of 
treatments might look like this: 
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(1) 
a 
e 


de 
f 
df 
ef 
def 


The alias subgroup of this 2°-* includes AD, BE and CF. If the assumptions 
that A, B and C do not influence y, and that D, EF and F do not influence y, 
are all correct, there is no confounding here. But the results will not give any 
clue to the failure of one or more of these assumptions. The two factorials can 
be adjoined in many ways. It will be better to look for a combination that is 
less sensitive to these assumptions. 

Since for eight runs the combined experiment must be a 2°~*, there will be 
seven interactions in the alias subgroup and there will be twenty-four letters. 
The average length of the interactions is therefore 3 3/7 letters, and so some 
interactions must be three-letter ones. There must then be some confounding 
of main effects with two-factor interactions. The alias subgroup 


-ADE ABDF —ACF 
—BEF  ABCE I 
—BCD CDEF 


guarantees that main effects on y, are confounded only with two-factor inter- 
actions on Y2, or on y, and y, together. (The letters D, EH, F are underlined.) 
Also estimates of main effects of factors D, FE and F on y, are only confounded 
with two-factor interactions on y,, or with mixed two-factor interactions. 

A simpler picture of what confounding is taking place is given by listing the 
sums of effects and 2 f.i. (two-factor interactions) that will be estimated by 
each of the familiar orthogonal contrasts: 


A — DE — CF D— AE — BC AB + CE + DF 
B-—EF-—-CD E-—-AD-—BF 
C—-BD-— AF F-—BE —AC 


The seven contrasts in the y, measurements will be given the interpretations 
A, B, C, BC, error, AC and AB respectively. The same contrasts in the y, 
measurements are given the interpretations DZ, EF, error, D, LE, F and DF 
respectively. The eight corresponding treatment combinations or experimental 
conditions are (1), aef, bde, abdf, acd, cdef, abce, bef, generated here by the three 
underlined combinations. 

The influence pattern 2-1-2 (meaning that two factors influence y, only, 
that one factor influences both responses, and that two factors influence y, 





‘tions 
he Y2 
1 DF 
ental 
three 


only, 
Ice Yo 


PARALLEL FRACTIONAL REPLICATES 265 


only) might also be studied in eight runs in a 2°’ with alias subgroup — ACE, 
—BCD, ABDE, I. The confounding pattern is given more explicitly by: 


Yo 
A —CE —CE 
B—CD —CD 
C — AE* — BD* C 
D — BC D 
E— AC E 
AB + DE DE 
AD* + BE* error error. 
The starred interactions are assumed to be zero. The eight experimental condi- 


tions are (1), abc, ae, bce, bd, acd, abde, cde, The seven contrasts among the 


y, measurements will be given the interpretation shown in the second column 
above; the same contrasts in the y, measurements will be given the names 
shown in the third column. 

For influence pattern 1-3-1 a plan in eight runs is possible although the 
2 f.i. will be confounded in pairs for each response. The alias subgroup J, ABCD, 
BCDE, AE produces the experimental conditions (1), abe, bc, ace, cd, abcde, 


bd, ade. The seven sums of confounded effects are now: 
A+E AB + CD + BE 
B + error AC +BD+CE 
C + error AD + BC + DE. 
D + error 


Although these plans are still too small for everyday use, Table 1 shows 
three two-factor-interaction-clear (2 f.i.c.) plans and four that confound 2 f.i. 
in pairs. These plans are not robust against wrongly assumed influence patterns 
and they give little or no information about error. 

Table 2 summarizes some 16-run, symmetrical influence pattern, bi-variate, 
2 f.i.c. plans. The column headed 3-1-3 is expanded here to show in more detail 
what confounding occurs and what sets of experimental conditions are required. 
It is known in this case that three factors, A, B and C affect only y; , that one 
factor, D, affects both y, and y, , and that three other factors, F, F and G affect 
only y, . Multiplying out all products of the three generators of the alias subgroup 
gives 


1. ABCG 1. X 2. CDFG 1. X 2. X 3. AEFG 
2. ABDF 1. X 3. BDEG 
3. ACDE 2. X 3. BCEF 
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TABLE 1 


Eight-run, symmetric-influence-pattern, bivariate, 2?-¢ plans 
I. Two-factor interaction clear 
II. Main effects clear, 2 f.i. confounded in pairs 


I. 
Influence 
Patterns 3-0-3 


Factor Letters ABC—DEF 


2 6 
q 3 
Members of (ABDF) . 
Alias Subgroups (ABCE) * 
(CDEF) 
(—ACF) 
Generators of acd 
Principal aef 
Blocks bef 
II, 
Influence 
Patterns 2-2-2 3-1-3 4-0-4 


Factor Letters AB-CD-EF ABC-D-EFG ABCD—EFGH 
P 6 7 8 
q 3 4 5 
Members of AE AE (ABCD) 
Alias Subgroups _ BF BF (EFGH) 
oa ca 
ABCD DH 
(DEFG) ABGH 
Generators of abcd ade abef 


Principal ace abef bdf aceg 
Blocks ade acf cdg adeh 


* Generators of alias subgroups underlined; other important members in parentheses. 


Thus, since every member of the alias subgroup has four letters, each main 
effect is confounded only with three-factor interactions. There will be three 
degrees of freedom for error estimation in each set of response contrasts, namely, 
those measuring D, E and F in y, and those measuring A, B and C in y, . The 
twenty-one 2 f.i. are confounded in the seven strings: 


AB + CG* + DF AF* + BD + EG 
AC + BG* + DE AG* + BC + EF 
AD+ BF*+CE* BE* + CF* + DG. 
AE*+CD + FG 
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TABLE 2 
Sixteen-run, symmetric-influence-pattern, two-factor-interaction-clear, bivariate, 29-2 plans 


Influence Patterns 1-3-1 2-2-2 3-1-3 4-0-4 


Factor Letters A-BCD-E AB-CD-EF ABC-D-EFG ABCD--EFGH 
p 5 6 7 8 
q 1 2 3 4 


Generators of —ABCDE ABCF ABCE 
Alias Subgroups ADEF ABDF 
ACDG 
BCDH 


Generators of abcf abcg abce 
Principal Blocks ace abdf abdf 
adef acde acdg 
af adg bedh 


The nine starred interactions are assumed negligible since their two factors 
taken together are stated not to influence y, or y2 . 

The set of 16 experimental conditions are produced from the four “generators 
of the principal block” by writing down all their products: 


l.abeg 1.X2.cdfg 2. X 4. bfg 1.3. xX 4. abe 
2.abdf 1.X3.bdeg 3. X 4. ceg 2.X3. xX 4. abcdefg 
3.acde 1.X4.bed 1.X2.X 3. aefg 1.2. X 3. X 4. def 
4.adg 2.X3.beef 1.X 2. X 4. acf (1). 
These combinations are easier to see when written in alphabetical order, thus: 
abcdefg bed 
abeg beef 
abdf bdeg 
abe bfg 
acde 
acf ceg 
adg def 
aefg (1). 


This plan is still not very robust against a wrongly assumed influence pattern. 
There are only four “error” degrees of freedom for each property, whose ir- 
regularity might indicate such a failure. However the improvement over the 
eight-run plans is noticeable. Clearly 32-run plans will give more “error” degrees 
of freedom, some of which may show significant inflation in case of serious 
mistakes in the assumed influence pattern. 
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Unsymmetrical influence patterns are common enough to be worth study. 
If the pattern is very skew, e.g. 4-1-1, then the requirements for y, , the re- 
sponse with the largest number of factors, dominate the plan. But moderate 
asymmetry can sometimes be adjusted nicely to the corresponding asymmetry 
in confounding of many 2?“ plans. 

Extensions to three or more properties do not seem easy to classify. They 
must be worked out in detail for each case. The simplest example is our old 
friend the 2°’, used now under the (unlikely) set of assumptions that only A 
and B affect y, , that only B and C affect y. , and that only A and C affect y; . 
The 2°-* may be viewed as three complete 2” factorials, one in A and B, one 
in B and C, one in A and C. 

A convenient way is needed to summarize the experimenter’s knowledge and 
feelings about the effects of each of K factors on each of R responses. A KX R 
“influence matrix” has been useful both in aiding the statistician to understand 
the limitations and advantages of the experimenter’s technical background, 
and to record the experimenter’s state of belief before the new round of experi- 
ments is started. Entries of —1, 0, +1 can be used to indicate the experimenter’s 
opinions about the sign and magnitudes of real effects. An 7 can be used to 
indicate ignorance. 

Useful proposals for the analysis of multivariate experiments seem to me 
to wait upon a more explicit classification of the types of responses to be measured 
in each run. Thus some responses are of similar properties, measured in the 
same units. Tensile strength before and after aging might be an example of 
such bi-variates. The second type of paired responses are similar or at least 
related physical or chemical properties, not of the same dimensions. Flexural 
strength and flexural modulus are a pair of this kind. Third might come proper- 
ties that do not depend on the same basic attributes of the material being 
produced. Brittleness and electrical conductivity, flavor and color, corrosion 
resistance and machinability are pairs that come in this category. We must 
expect different sorts of dependence, functional and probabilistic, for these 
classes of responses. A single statistical model will hardly apply to all. 

Detailed numerical examples will be published later. 
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TECHNOMETRICS May, 1960 


NOTICES 


DEPARTMENT OF STATISTICS 
Tue Fioripa State UNIVERSITY 


To assist in meeting the demands for well trained statisticians, the Florida 
State University in 1959 founded a new Department of Statistics with immediate 
authorization to offer programs leading to the Bachelor of Science and Master 
of Science degrees in statistics. With the further development of the department, 
programs leading tothe Doctor of Philosophy degree are anticipated. 

The Bachelor of Science degree in statistics is designed to prepare the student 
for continuing graduate work, for employment related to statistical compu- 
tations, and for employment with statistics groups in government, industry 
and business. A major in statistics consists of a minimum of eighteen semester 
hours in statistics courses at the junior and senior college level. Competency 
in mathematics is essential and a minimum of eighteen semester hours in mathe- 
matics is required. Survey courses in areas of application of statistics are recom- 
mended and at least twelve semester hours are to be completed in such an area. 

The Department of Statistics offers the Master of Science degree in statistics 
with a thesis requirement. Prior work in statistics is not a requirement for 
admission to graduate work, but mathematical training comparable to an under- 
graduate major in mathematics is necessary. Graduates in the sciences may 
be admitted to graduate work in statistics, but it may be necessary to make 
up deficiencies in mathematics. The program for the Master of Science degree 
has been developed so that students completing the program may continue 
directly in graduate training in statistics toward the Doctor’s degree at other 
universities. The Department of Statistics at the Florida State University offers 
unusual individual attention and cooperation with the beginning graduate 
student. 

In addition to degree programs, the Department of Statistics is engaged 
in teaching and consultation in statistics with many other university depart- 
ments. Courses are offered for other curricula. Faculty and students consult 
on research projects and theses. Members of the Department of Statistics are 
engaged in research in statistics in a variety of areas of interest. The editorial 
office of Bromerrics, the journal of the International Biometric Society is 
in the Department of Statistics. 

For complete details on degree requirements, the catalog bulletins of the 
University should be consulted. 

ASSISTANTSHIPS. Graduate fellowships are available at stipends for the 
academic year of $1500 for a Master’s Degree candidate. These awards require 
no service and provide a full academic year of residence. Fellows pay the regular 
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resident registration fee, but out-or-state tuition is waived. Applications for 
fellowships should be filed by February 15th for the following academic year. 
Later applications will be considered if funds are still available. Fellowships 
are tax exempt. 

Graduate assistantships provide stipends for the academic year varying from 
$160 to $240 per month for half-time service. Both teaching and research 
assistantships are available in the Department of Statistics and in association 
with the Computing Center and other departments. The amount of the stipend 
is determined by the qualifications of the candidate and the nature and the 
amount of the part-time services to be rendered. Assistants pay the regular 
resident registration fee, but out-of-state tuition is waived. Applications for 
assistantships are considered at the time received, but should normally be 
filed by March 15 for the following academic year. 

Additional information on fellowships and assistantships may be obtained 
from the Department of Statistics and from the Graduate School. Applications 
for fellowships are made to the Graduate School and for assistantships to the 
Department of Statistics. 


COURSES. Courses in statistics offered by the Department of Statistics 
are as follows: 


For Undergraduates: 


216 Elementary Statistics 
316 General Theory of Statistics 


For Advanced Undergraduates and Graduates: 


406 Statistical Procedures for Research 
407 Statistical Methods and Design 
426 Fundamentals of Sampling 

431 Mathematical Statistics I 

432 Mathematical Statistics II 

446 Theory of Probability 

491-2-3 Directed Individual Study 


For Graduates: 


501 Theory of Statistics I 

502 Theory of Statistics IT 

506 Advanced Topics in Statistical Methods 
507 Nonparametric and Rank Order Statistics 
508 Introduction to Multivariate Analysis 
509 Advanced Statistical Inference 
591-A*B*C Directed Individual Study 

599 Master’s Thesis 


Descriptions and prerequisites for courses with other details are given in the 
Florida State University general catalogue and Graduate Bulletin. 
Additional information and assistance with applications may be obtained 
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by writing to Dr. Ralph A. Bradley, Head, Department of Statistics, Florida 
State University, Tallahassee, Florida. Telephone, Tallahassee 2-8000, Exten- 
sion 2369. 





CorNELL UNIVERSITY INDUSTRIAL 
ENGINEERING SEMINARS 


The Department of Industrial and Engineering Administration, Sibley School 
of Mechanical Engineering, Cornell University, announces the ‘Cornell Uni- 
versity Industrial Engineering Seminars” to be held at Cornell University, 
Ithaca, New York, June 14 through June 17, 1960. Participants may enroll 
in one of the following seven groups and also attend general sessions. 


Industrial Management 

Engineering Administration 

Operations Management of the Smaller Company 

Work Measurement 

Systems Simulation Using Digital Computers 

Statistical Decision-Making: Theory and Applications 
Statistical Reliability Analysis: Theory and Applications 


Specialists from both industry and the staff of Cornell University will provide 
the seminars. The session is planned for operating management personnel in 
line supervision and staff positions in industrial engineering, production engineer- 
ing, engineering administration, research and development, quality control, 
production control, cost control, materials control, purchasing, and data pro- 
cessing. 

For additional information write J. W. Gavett, Seminars Coordinator, Depart- 
ment of Industrial and Engineering Administration, Upson Hall, Cornell Uni- 
versity, Ithaca, New York. 


STATISTICAL QuALITY CONTROL FOR THE TEXTILE INDUSTRY 


The School of Textiles of the North Carolina State College, Raleigh, North 
Carolina, in cooperation with the College Extension Division will sponsor this 
short course May 2-13, 1960. 


CoursE CoNnTENT 
First Week, May 2-6, 1960 
Monday 
Registration and Introduction 
Populations and Samples 
Frequency Distributions 


Methods of Calculation and an of Standard Deviation 
Practice Problems 





272 NOTICES 


Tuesday 


Calculations of Coefficient of Variation and Applications to Textile Processes 
Statistical Quality Control Charts (X and R) 
Practice Problems and Case Studies 


Wednesday 


Statistical Quality Control Charts (P, PN and @) 
Applications, Case Studies and Practice Problems 


Thursday 

Quality Control Systems by Areas for Textile Plants 

Interpretation and Application of Testing Data to Manufacturing Processes 
Friday 

Applications of Quality Control Techniques with Case Studies 

Review with Question and Answer Period 


Second Week, May 9-13, 1960 


Monday 


Tests for Significance of Differences and Confidence Intervals Applied to 
Textiles 


Tuesday 


Continuation of the Application of Significance Testing and Confidence 
Intervals to Textiles, with Practical Problems 
Correlation and other Relationships of Variables as Related to Textiles 


Wednesday 


Applications of Correlation to Textile Problem 
Analysis of Variance 


Thursday 
Case Studies in the Use of Analysis of Variance 
Sampling Theory as Applied to Continuous Production and Lot-by-Lot 
Inspection 
Continuation of Sampling Theory and Special Applications of Statistical 
Techniques to Textile Quality Control Problems 
The course will begin on Monday, May 2 and end on Friday, May 13. There 
will be seven hours of instruction per day which will include lectures, case 
studies, and practical problems. 


In order to provide personal instruction, enrollment is limited to twenty. 
Registration for the course will be accepted in the order of paid registrations 
received. 

Each participant for the full course will be provided with a notebook con- 
taining the necessary materials, as well as a copy of the appropriate textbook. 
Each enrollee is encouraged to have a slide rule for use in the course. 
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The instruction for the course will be under the direction of Prof. D. S. Hamby 
of the Department of Textile Technology, School of Textiles, N. C. State College. 
Assisting as full-time instructors are Dr. R. J. Hader, Department of Experi- 
mental Statistics, N. C. State College, and Prof. W. C. Stuckey, Jr., Department 
of Textile Technology, School of Textiles, N. C. State College. In addition, 
guest lecturers will be invited for discussions during the case studies. All of 
the instructors have had experience in the teaching of statistical quality control, 
both to students and personnel in short courses. They are also well qualified 
in industrial applications because of their work as consultants to numerous 
firms in the textile industry. 

Registration will be conducted in the Nelson Building (School of Textiles) 
from 9:00 to 9:30 A.M. on Monday, May 2, 1960. The registration fee is $160.00 
per person for the full course. This fee encludes tuition, supplies, and text 
materials, but does not include housing accommodations or meals. 


Tue SraTisTICAL SUMMER SEMINAR 


The Statistical Summer Seminar will hold its tenth meeting July 18-25, 1960 
at the NUY Gould House, Ardsley-on-Hudson, New York to discuss ‘The 
Use of Computers in the Sixties—Boom or Bust?’ Dr. Meyer A. Efroymson 
of Esso Research and Prof. Fred C. Leone of Case Institute of Technology are 
organizing a program which will include several noted speakers. 

This seminar, a private non-profit organization, was organized in 1948 by 
Profs. Frederick Mosteller and John Tukey to explore new fields of application 
for statistics. At present it is under the direction of Dr. Milton E. Terry of 
the Bell Telephone Laboratories. 

The registration fee of $175.00 includes comfortable living accommodations. 
Grants in aid are available to government and academic personnel. For further 
information write to Dr. Ramanathan Gnanadesikan, Bell Telephone Labora- 
tories, Murray Hill, New Jersey. 


Puup AND Paper Sratistics Cours—E ANNOUNCED 


The ninth Pulp and Paper Statistics Course will be held August 8th through 
19th, 1960, at Brevard College, Brevard, N. C., and will be jointly sponsored 
by the Institute of Statistics of North Carolina State College, The Technical 
Section of the Canadian Pulp and Paper Association, and the Technical Associ- 
ation of the Pulp and Paper Industry (TAPPI). 

Two courses of instruction, one basic and one advanced, will be held con- 
currently during the two-week period. The courses will be concerned with the 
techniques of statistical analysis in research, control and management methods 
in the pulp and paper industry. 

The Elementary Course is intended for those people who wish to acquire 
some knowledge of, and facility with, the simpler but exceedingly valuable 
statistical techniques and who have had little or no previous formal training 
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in statistical methods. The training received from this course should prove 
profitable to organizations now making little or no use of statistical methods 
but wishing to know more about how they work and how to use them. It will 
also be profitable to organizations with statistical programs already operating 
who wish to strengthen their position by having additional men with statistical 
training available within their organization. 

The Advanced Seminar is intended for people who have experience and definite 
facility in the use of statistical techniques and who wish to gain training and 
practice in special problems that arise in application. The course will consist 
of laboratory sessions devoted to the consideration of specific problems in the 
industry, and the practical analysis of these problems. 

Further information and application forms for the course may be obtained 
from the Institute of Statistics, North Carolina State College, Raleigh, N. C. 





THE ANNALS OF 
MATHEMATICS STATISTICS 


Vol. 31, No. 1—March, 1960 


Contents 
The Algebra of a Linear Hypothesis Henry B. Mann 
An Operational Approach to the r-Way Crossed Classification ..J. D. Bankier 
Second Order Rotatable Designs in Four or More Dimensions 
Norman R. Draper 
A Matrix Substitution Method of Constructing Partially Balanced Designs 
B. V. Shah 
Some Approximations to the Binomial Distribution Function ..R. R. Bahadur 
On the Mixture of Distributions Henry Teicher 
On Interchanging Limits and Integrals John W. Pratt 


Moments of the Absolute Difference and Absolute Deviation of Discrete 
ENN NE) acs dy: dil aia ase pee Rn wad & rata eee sasieloues S. K. Katti 


Priority Queues Rupert G. Miller, Jr. 
On Time Dependent Queuing Processes ........ J. Keilson and A. Kooharian 


A Geometry of Binary Sequences Associated with Group Alphabets in 
Information Theory R. C. Bose and Roy R. Kuebler, Jr. 


A Necessary and Sufficient Condition for the Existence of Consistent 
Estimates Lucien LeCam and Lorraine Schwartz 


The Distribution of the Latent Roots of the Covariance Matrix . Alan T. James 
Two-Sample Tests for Multivariate Distributions Lionel Weiss 


A Modification of the Sequential Probability Ratio Test Analysis to Reduce . 
the Sample Size T. W. Anderson 
Sequential Tolerance Regions Sam C. Saunders 


A Useful Generalization of the Stein Two-Sample Procedure .R. Wormleighton 


Notes: 


Sums of Small Powers of Independent Random Variables ..J. M. Shapiro 
On the Median of the Distribution of Exceedances K. Sarkadi 


The purpose of the Institute of Mathematical Statistics is to encourage the development, 
dissemination and application of mathematical statistics. Membership dues, which include a 
subscription to the Annals of Mathematical Statistics are $10.00 per year for residents of the 
United States of America and $5.00 a year for residents of other countries. Inquiries regarding 
membership to the Institute should be sent to the Secretary: G. E. Nicholson Jr., Department 
of Statistics, University of North Carolina, Chapel Hill, North Carolina, USA. 
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CONTENTS 
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BARTLETT, M. S., Gower, J. C. & Lesiiz, P. H. A comparison of theoretical and 
empirical results for some stochastic population models 
KENDALL, Davin G. Birth-and-death processes, and the theory of carcinogenesis 


AsHForD, J. R., Smit, C. S. & BRown, SUSANNAH. The quantal response analy- 
sis of a series of biological assays on the same subject 


WALKER, A. M. Some consequences of superimposed error in time series analysis 


Fincu, P. D. Deterministic customer impatience in the queueing system 
GI/M/1 


Barton, D. E., Davin, F. N. & Fix, Evetyn. The polykays of the natural 
numbers 


Srivastava, A. B. L. The distribution of regression coefficients in samples from 
bivariate non-normal populations I. Theoretical investigation 


KASTENBAUM, Marvin A. The separation of molecular compounds by counter- 
current dialysis: a stochastic process 


Saw, J. G. A note on the error after a number of terms of the David-Johnson 
series for the expected values of normal order statistics 


Watson, G. S. More significance tests on the sphere 


Jounson, N. L. An approximation to the multinomial distribution: more 
properties and applications 


VAGHOLKAR, M. K. & WETHERILL, G. B. The most economical binomial sequen- 
tial probability ratio test 


McGrecor, J. R. An approximate test for serial correlation in polynomial 
regression 


BERKSON, JOSEPH. Nomograms for fitting the logistic function by maximum 
likelihood 


HAIGHT, FRANK A. & BREUER, MELVIN ALLEN. The Borel-Tanner distribution 
Burr, E. J. The distribution of Kendall’s score S for a pair of tied rankings 


Miscellanea: 


Contributions by D. H. Buate, C. R. BLytH & C. L. Curme, C. Burrows, C. W. CLUNIES- 
Ross & R. H. RiIFFENBURGH, P. D. Fincu, R. G. LAHA & E. Lukacs, P. H. Leste, H. D. 
PATTERSON, B. R. Rao, G. L. Turin, C. A. Wirkins, M. E. WIsE. 


Reviews: 

Other Books received: 
The subscription, payable in advance, is 54/- (or $8.00) per volume (including postage). 
Cheques should be made payable to Biometrika, crossed “A/c Biometrika Trust” and 


sent to the Secretary, Biometrika Office, University College, London, W. C. 1. All foreign 
cheques must be drawn on a Bank having a London agency. 
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Vol. 16, No. 2 CONTENTS June 1960 


The Role of Statistics in Physiological Research C. W. Emmens 


Rank Analysis of Incomplete Block Designs: Method of Paired 
Comparisons Employing Unequal Repetitions on Pairs Otto Dykstra 


The Treatment of Reciprocal Interaction, with or without Lag, 
in Path Analysis Sewall Wright 


Estimating the Parameter in a Conditional Poisson 
Distribution A. C. Cohen, Jr. 
Models for the Interpretation of Experiments Using Tracer 
Compounds Jerome Cornfield, 
Jesse Steinfield, and 
Samuel W. Greenhouse 
Genetic Correlations with Multiple Alleles R. G. Stanton 


Analysis of Experiments Measuring Threshold Taste E. K. Harris 


An Unbiased Sampling and Estimation Procedure for Creel 
Censuses of Fishermen D. 8S. Robson 


The Mathematical Theory of Biological Assay of a Local 
Anaesthetic N. K. Chakravarti 


t 
QUERIES ann NOTES 
Calculation of Inbreeding in Family Selection Studies on the 
IBM 650 Data Processing Machine K. Hoen and 


A. H. E. Grandage 
A Note on Optimum Family Size in Selection Programs Allan Robertson 
On Heuristic Estimation Methods Malcolm E. Turner 


A Genetic Application of the Schumann-Bradley Table Stanley Wearden 


Biometrics is published quarterly by The Biometric Society. Its objects are to describe 
and exemplify the use of mathematical and statistical methods in biological and related 
sciences, in a form assimilable by experimenters. The annual non-member subscription rate 
is $7. Inquiries, orders for back issues and non-member subscriptions should be addressed to: 
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PREPARATION OF MANUSCRIPTS 


Manuscripts should be submitted to the office of the editor: J. S. Hunter, 
Mathematics Research Center, U. 8. Army; The University of Wisconsin; 
Madison, Wisconsin. Each manuscript should be typewritten, double spaced, 
with wide margins at sides, top, and bottom. The original should be submitted 
with two additional copies, on paper that will take corrections. Dittoed or 
mimeographed papers are acceptable only if completely legible. Footnotes 
should be avoided and replaced by remarks in the text, or placed in an appendix. 
Preferably, references in the manuscript should appear as (Jones, A. B., 1958), 
and again later in alphabetical order in a list of references. Alternatively, refer- 
ences may be numbered, e.g. [1], as they appear in the manuscript and be listed 
in this sequence in the list of references. In the reference list, each reference 
should contain, in the order indicated, the name and initials of the author 
followed by those of the co-authors, date of publication, title of reference, 
souree, volume number and page. References to books should include pub- 
lisher’s name and location. 

Figures, charts, and diagrams should be professionally drawn on plain white 
paper or tracing cloth in black India ink twice the size they are to be printed. 
A full page diagram, in print, measures 7.25 4.75 inches. 

As far as possible, formulas should be typewritten and symbols not available 
on a typewriter carefully inserted in ink. Authors are asked to keep in mind the 
typographical difficulties of complicated mathematical formulae. The difference 
between capital and lower-case letters should be clearly shown; care should be 
taken to avoid confusion between such pairs as zero and the letter O, the numeral 
1 and the letter J, numeral 1 used as superscript and prime (’), alpha and a, kappa 
and k, mu and u, nu and », eta and n, ete. Subscripts or superscripts should be 
clearly below or above the line. Bars above groups of letters (e.g., log x) and 
underlined letters (e.g., z) are difficult to print and should be avoided. Symbols 
are automatically italicized by the printer and should not be underlined on 
manuscripts. Boldface letters may be indicated by underlining with a wavy line 
on the manuscript; boldface subscripts and superscripts are not available. 
Complicated exponentials should be —. with the symbol exp particu- 
larly when appearing in the text, that is, 


exp [(a” + b”)'””] should be used in place of ecartente 


In writing square roots the fractional exponent is preferable to the radical sign. 
Fractions in the body of «he text (and when possible in displayed expressions) 
and fractions occurring in the numerators or denominators of fractions are 
preferably written with the solidus; thus 


a+b 
c+a 
Authors will ordinarily receive only galley proofs. Fifty offprints without 


covers will be furnished free. Costs for additional reprints and covers can be 
furnished on request. 
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